How old is our universe? Is the expansion of space accelerating? Or will it eventually collapse into a singularity? The evolution of our universe can, surprisingly enough, be modeled by solving a set of quite simple ordinary differential equations (ODEs) called the Friedmann equations.

When viewed on a large enough scale, our universe seem to be both homogeneous and isotropic. It is uniform in distribution of matter and radiation in all directions, and has no observable irregularities. Together with the assumption that the spacetime evolves as a perfect fluid, Alexander Friedmann derived his set of equation from Einstein's field equations of gravity in 1922. The equations depends on the density and pressure of the matter, radiation, and vacuum components of the universe.

The variable $a(t)$ is called the *scalar factor* and is a parameterization of the relative expansion of the universe. It is time dependent and its current value ($t = 0$) is defined as $a(0) = a_0 = 1$. Subscripts '0' denote the quantities evaluated today. $a$ increases with the expansion of the universe. Thus, at the Big Bang singularity, $a(t_{BB}) = 0$ with $t_{BB} < 0$.

The first Friedmann equation is a function of the scalar factors time derivative $da/dt = \dot{a}$,

\begin{equation} \Big(\frac{\dot{a}}{a}\Big)^{2} = \frac{8\pi G}{3}\rho - \frac{\kappa}{a}^2 . \label{first} \end{equation}

The second of its acceleration:

\begin{equation} \frac{\ddot{a}}{a} = -\frac{4\pi G}{3}(\rho + 3P) . \label{second} \end{equation}

Here $G$ is the gravitational constant. If you want to learn how the equations are derived from Einstein's field equations of gravity, see [1].

The parameter $\rho = \sum_i \rho_i$ is the energy density of the perfect fluid. The different components $\rho_i$ have different contributions to the pressure, $P$, given by the equation of state $w = P/\rho$, $$P = \sum_i P_i = \sum_i w_i\rho_i. $$

The different components are:

\begin{equation} \rho \propto a^{-3(1 + w)} = \begin{cases} a^{-3} & \mathrm{matter} \\ a^{-4} & \mathrm{radiation} \\ a^{0} & \mathrm{vacuum} \end{cases} \label{rho} \end{equation}

The vacuum density corresponds to the mysterious *dark energy* which exerts negative pressure on the universe. As $a$ increases, the value of this component will eventually dominate the expansion.

The parameter $\kappa$ denotes the curvature of the isotropic and homogeneous three dimensional space. In terms of two dimensions, positive curvature ($\kappa = 1$) has spherical spatial geometry, negative curvature ($\kappa = -1$) has hyperbolical spatial geometry, and no curvature ($\kappa = 0$) has flat spatial geometry.

The first Friedmann equations can be rewritten in terms if the Hubble parameter

$$ H \equiv \frac{1}{a}\frac{da}{dt} = \frac{\dot{a}}{a}, $$

and the dimensionless density parameters

$$ \Omega_{I, 0} \equiv \frac{\rho_{I, 0}}{\rho_{crit, 0}} ,$$

where the critical density

$$\rho_{crit, 0} \equiv \frac{3 H_0^2}{8\pi G} ,$$

and the lower index $I \in \{m, r, \Lambda\}$ denotes the specific density component, corresponding to matter, radiation and vacuum energy, respectively. $H_0$ denotes the Hubble parameter today, often referred to as the Hubble constant.

Using these parameter with the relations \eqref{rho}, Friedmanns first equation becomes

\begin{equation} H^2(a) = \Big(\frac{\dot{a}}{a}\Big)^2 = H_0^2 \Big{[} \Omega_{r, 0}\Big(\frac{a_0}{a}\Big)^{4} + \Omega_{m, 0} \Big(\frac{a_0}{a}\Big)^{3} + \Omega_{\kappa, 0} \Big(\frac{a_0}{a}\Big)^{2} + \Omega_{\Lambda, 0}\Big{]} , \label{firstRewritten} \end{equation}

where the curvature density parameter is defined as $\Omega_{\Lambda, 0} \equiv -\kappa/(a_0 H_0)^2$.

From the 2018 results of the observational spacecraft Planck [2] and [1], the cosmological parameters of our universe were presented to be:

\begin{equation} \begin{split} \Omega_m &= 0.315 \pm 0.007 \\ \Omega_\Lambda &= 0.6847 \pm 0.0073 \\ \Omega_r &= 9.4 \cdot 10^{-5} \\ \mid \Omega_\kappa \mid &\leq 0.01 \\ H_0 [\mathrm{km}\text{ }\mathrm{s}^{-1} \mathrm{Mpc}^{-1}] &= 67.4 \pm 0.5 . \end{split} \end{equation}

In [1]:

```
omega_m = 0.315
omega_lambda = 0.6847
omega_r = 9.4e-5
omega_kappa = 1 - omega_m - omega_lambda - omega_r
H0 = 67.4 # [Km/(s Mpc)]
H0 = H0*1000*(3600*24*365*1e9)/3e22 # [1/Gyr] (Gyr = 10^9 years)
```

To estimate the age of the universe, the equation \eqref{firstRewritten} can be rearranged to an integral

\begin{equation} H_0dt = \int_0^{a_0} \frac{da}{\Big{[} \Omega_{r, 0}a^{-2} + \Omega_{m, 0}a^{-1} + \Omega_{\kappa, 0} + \Omega_{\Lambda, 0}a^2\Big{]}^{1/2}}, \label{age} \end{equation}

where $a_0$ = 1.

This integral can be solved using the function `scipy.integrate.quad`

.

In [2]:

```
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import quad
import matplotlib.pyplot as plt
import sys
%matplotlib inline
```

In [3]:

```
# Defining the integrand as a function
def integrand(a, omega_r, omega_m, omega_kappa, omega_lambda):
return (omega_r*a**(-2) + omega_m/a + omega_kappa + omega_lambda*a**2)**(-1/2)
# Solve definite integral with limits (0, 1) by using scipy.integrate.quad
H0t = quad(integrand, 0, 1, args=(omega_r, omega_m, omega_kappa, omega_lambda))
t = H0t[0]/H0
print("Age of universe [Gyr]:\t%.1f"%(t))
```

This result is fairly close to the current measurement by Planck of about 13.8 billion years.

Combining \eqref{first} with \eqref{second} yields the system of ODE's

\begin{equation} \begin{cases} & \dot{a} = \frac{da}{dt} \\ &\frac{d\dot{a}}{dt} = -\frac{H_0^2}{2}\Big(2\Omega_{r, 0}a^{-4} + \Omega_{m, 0}a^{-3} - 2\Omega_{\Lambda, 0}\Big)\Omega_{m, 0}\Big(\frac{\dot{a}^2}{H_0^2} - \Omega_{r, 0}a^{-2} - \Omega_{\Lambda, 0}a^2 - \Omega_{\kappa, 0}a_0^2\Big)^{-1} \label{dderivative} \end{cases} \end{equation}

With initial conditions

\begin{equation} \begin{split} &a_0 = 1 \\ &\dot{a}_0 = H_0\sqrt{\Omega_{m, 0}a_0^{-1} + \Omega{r, 0}a_0^{-2} + \Omega{\Lambda, 0}a_0^2 + \Omega_{\kappa, 0}a_0^2} . \end{split} \end{equation}

Let's define a function for calculating the derivatives of $a$.

In [4]:

```
def friedmann(y, t, omega_r, omega_m, omega_lambda, omega_kappa, H0, a0):
"""
Returns the first and second derivative of the scalar factor,
given the the current value of the scalar factor and its first
derivative.
Parameters:
y 1D array with two elements
t 1D array (necessary for .odeint function)
Omega_x Float. Density parameters
H0 Float. Hubble parameter
a0 Float/integer. Scalar factor today
Returns:
dydt 1D array with two elements
"""
# Defining non- and first derivative based on input-
a, adot = y
dydt = [adot, -H0**2/2*(omega_m*a**(-3) + 2*omega_r*a**(-4) -
2*omega_lambda)*omega_m*(adot**2/H0**2 - omega_r*a**(-2) -
omega_lambda*a**2 - omega_kappa*a0**2)**(-1)]
return dydt
```

Next, we solve the system of equations using the integrator `odeint`

from the `scipy.integrate`

library. We choose four different sets of density parameters to get different evolution of the scalar factor. The values are not chosen randomly, but in order to achieve specific models corresponding to theoretical scenarios. The system is integrated both backwards and forwards to get the evolution both in the past and the future.

In [7]:

```
# Defining constant parameters
H0 = 100 # Hubble constant [Km/(s Mpc)]
H0 = H0*1000*(3600*24*365*1e9)/3e22 # Hubble constant [1/Gyr] (Gyr = 10^9 years)
a0 = 1 # Scalar factor today
# Limits of integration (billion years)
up_lim = 30
low_lim = -30
# Resolution of integration
res = int(1e4)
# Sequence of time steps
t_future = np.linspace(0, up_lim, res)
t_past = np.linspace(0, low_lim, res)
t = np.concatenate([np.flip(t_past, 0), t_future[1:]])
# Defining four different sets of density parameters values
# omega = [omega_r, omega_m, omega_lambda, omega_kappa]
omega1 = [1e-4, 0.3, 0, 1-1e-4-0.3-0]
omega2 = [1e-4, 0.3, 0.7, 1-1e-4-0.3-0.7]
omega3 = [1e-4, 5, 0, 1-1e-4-5-0]
omega4 = [1e-4, 1, 0, 1-1e-4-1-0]
SETS = [omega1, omega2, omega3, omega4]
# Initializing matrix to hold solutions
evolution = np.zeros((len(SETS), 2*res-1))
for i in range(len(SETS)):
omega_r, omega_m, omega_lambda, omega_kappa = SETS[i]
# Initial value of the time derivative of the scalar factor
adot0 = H0*np.sqrt(omega_m/a0 + omega_r*a0**(-2) + omega_lambda*a0**2 + omega_kappa*a0**2)
# List of initial conditions
y0 = [a0, adot0]
# Find numerical solutions
sol_future = np.zeros((2, res))
sol_past = np.zeros((2, res))
sol_future = odeint(friedmann, y0, t_future, args=(omega_r, omega_m, omega_lambda, omega_kappa, H0, a0))
sol_past = odeint(friedmann, y0, t_past, args=(omega_r, omega_m, omega_lambda, omega_kappa, H0, a0))
# Merge the solution for past and future and insert in matrix
evolution[i, :] = np.concatenate([np.flip(sol_past[:, 0], 0), sol_future[:, 0][1:]])
# Plot the results
plt.figure(figsize=(12, 8), dpi=300)
for i in range(len(SETS)):
plt.plot(t, evolution[i],
label=(r'$\Omega_{r, 0} = %.4f, \Omega_{m, 0} = %.1f, \Omega_{\Lambda, 0} = %.1f,'
'\Omega_{\kappa, 0} = %.1f$'%(SETS[i][0], SETS[i][1], SETS[i][2], SETS[i][3])))
plt.xlabel("Time, Gyr")
plt.ylabel(r'Scalar factor, $a$')
plt.ylim(0, 4)
plt.xlim(-15, 30)
plt.title("Relative size of the universe as a function of time")
plt.legend(loc='best')
plt.grid(True)
plt.show()
```

The plot above shows the evolution of the scalar factor, and hence the relative size of the universe, for the different density parameters defined. All graphs meet at $a=1$, which is today. Therefore, our universe could have followed either of the four evolutionary paths modeled in the plot, depending on the true value of the parameters.

The green curve corresponds to a spherical universe with $\Omega_{\kappa, 0} < 0$ and thus the curvature $\kappa = 1$. Noticeably, it is the youngest universe, less than five billion years old. This universe will eventually recollapse into a singularity, a scenario referred to as the *Big Crunch*.

The red curve corresponds to the *Einstein- de Sitter Universe*. The curvature is flat and the universe expands eternally but is decelerating.

In the blue curve model, the curvature is hyperbolic ($\kappa = -1$). Also here the expansion is eternal but decelerating.

Finally, the orange. The curve which by current observations most resembles our own universe. The curvature density $\Omega_{\kappa, 0} \approx 0$, indicating a flat universe. Interestingly, the cosmological constant is significant, resulting in an accelerating expansion of the universe as $a > 1$.

[1]: Amsterdam Cosmology Group, *Lecture Notes: Chapter 1*, University of Amsterdam and Nikhef.

[2]: Planck Collaboration, *Planck 2018 Results. IV. Cosmological Parameters*, arXiv:1807.06209 [astro-ph.CO].

[3]: Cosmology, *Multi-Component Universes*, Department of Astronomy, University of Stockholm.