Most people are familiar with the standard model of cosmology: the universe (or at least the part of it we can see) began in a Big Bang and has been expanding and cooling off ever since. Mathematical models of the expansion of the universe have been around since the 1920s, but only in the last decade or so have measurements become precise enough to plug specific numbers into those models. For instance, just a few years ago we knew only enough to put the age of the universe between 10 and 20 billion years old, but today we know it is 13.75 billion years old, and that’s accurate to within 1%!

Models of the expanding universe are all based on general relativity. GR is conceptually not too complicated—“spacetime tells matter how to move, and matter tells spacetime how to curve”, as the saying goes—but the mathematical machinery required to make it all precise is famously arcane. Fortunately, models of the expanding universe turn out to actually be quite simple.

On the largest scales of space, the universe is homogeneous and isotropic. On these scales, galaxies are like molecules, and the matter in the universe can be well approximated as a continuous, uniform fluid filling all of space. This simplification brings GR from 10 partial differential equations defining 10 functions of space and time—a numerical nightmare—down to just one ordinary differential equation, defining one function of time, and this equation is easily integrated numerically.

This equation is called the Friedmann equation, and it is: $$\frac{d^2 a}{dt^2} = -\tfrac{4}{3}\pi G (\rho + 3p) a(t).$$ Here, $a(t)$ is the scale factor, a function of time that describes the size of the universe as a ratio to its size today. That is, if we take $t = 0$ to be the present day, then $a(0) = 1$. This is the function for which we’re trying to solve. The remaining values are parameters: $\rho$ (the Greek letter rho) is the energy density of the contents of the universe (its matter and energy); $p$ is the pressure of those contents, and $G$ is the gravitational constant.

We can’t solve this equation quite yet, because the density and pressure are actually functions of the scale factor. For example, as the universe grows larger, all matter gets spread further apart, so its density and pressure decrease. In order to solve the Friedmann equation, we need to pin down how $\rho$ and $p$ depend on $a$. To do this, we’ll need to separate the contents of the universe into different categories that have different density-vs-pressure curves.

There are three main kinds of stuff in the universe that need to be taken into account:

Matter. This category contains all stars, planets, black holes, nebulas, dust clouds, and dark matter—everything made of particles or objects that move slowly compared to lightspeed. (We don’t know what kind of particles dark matter is made of, but we’re pretty sure they’re nonrelativistic.) On a universal scale, matter has no significant pressure, so $p_M = 0$. It spreads out as the universe expands, inversely proportional to volume, so $\rho_M(t) = \rho_{M0} / [a(t)]^3$, where $\rho_{M0}$ is the present-day density of matter.

Radiation. In the context of cosmology, radiation refers to any particles that move at a relativistic (near-light) speed—mainly photons and neutrinos. Special relativity implies that these particles have a characteristic pressure of $p_R = \tfrac{1}{3} \rho_R$. Like matter, radiation also spreads out as the universe expands, but the particles it is made of also lose energy and momentum to the expansion, so $\rho_R(t) = \rho_{R0} / [a(t)]^4$—a stronger falloff than that of matter. Again, $\rho_{R0}$ is the present-day radiation energy density.

Dark energy (aka vacuum energy). We frankly don’t know what dark energy is made of or how it behaves, except that it has a negative pressure—that’s why it’s causing the universe to expand faster, instead of slowing down. Because of the uncertainty, physicists usually express its behavior in terms of an unknown parameter, $w$, where $p_{DE} = w\rho_{DE}$. This leads to a scale factor dependence of $\rho_{DE}(t) = \rho_{DE0} / [a(t)]^{3w+3}$. We know the value of $w$ is in the vicinity of −1, but to date this is only known to within about 10%. It is also not known whether $w$ is a constant, or whether it can change over time. The value of $w$ can make a big difference to the long-term fate of the universe, as we’ll see later.

Now that we have equations for the behavior of each component of the contents of the universe, we can put them back into the Friedmann equation, adding all the densities and pressures together: $$\frac{d^2 a}{dt^2} = -\tfrac{4}{3}\pi G \left[ \frac{\rho_{M0}}{a^2} + 2 \frac{\rho_{R0}}{a^3} + (1 + 3w) \frac{\rho_{DE0}}{a^{3w + 2}} \right].$$ That looks kind of complicated, but it’s now in a form that we can attack with numerical integration. We just need some data to start it off with.

As mentioned, we’ll use $t = 0$ as the present day, and since the scale factor is relative to the present, $a(0) = 1$. The remaining values at the present day are as follows (from here, with unit conversions where necessary):

• $da/dt = 7.20 \times 10^{-11} \;\text{yr}^{-1}$ (this is the Hubble “constant”—not really a constant on a universal time-scale, but constant enough from a human point of view!)
• $\rho_{M0} = 2.53 \times 10^{-27} \;\text{kg/m}^3$
• $\rho_{R0} = 5.6 \times 10^{-31} \;\text{kg/m}^3$
• $\rho_{DE0} = 6.78 \times 10^{-27} \;\text{kg/m}^3$

Notice that the densities of matter and dark energy at the present time are about the same order of magnitude, but the radiation density is 10,000 times lower. Radiation was very important in the early universe, but it hardly affects the expansion of the universe today, and going into the future it will be even more insignificant.

I used Python to integrate the Friedmann equation; I’ll link to the code at the end of this post. Starting from the present day and going both forward and backward in time, the scale factor looks like this (click to zoom):

There are several interesting things in this graph. First of all, notice that the scale factor goes to zero at −13.75 billion years. This precisely matches the standard value for the age of the universe computed by real physicists, which gives me confidence that this model is on the right track! Initially, the universe expanded quickly, but in the first couple billion years it slowed down—the attractive gravity of all the matter and radiation in it acted as a brake on the expansion. Then, sometime within the last five billion years, the expansion started to speed up again. This is happening because of dark energy, which has a repulsive gravitational effect due to its negative pressure. In the future, matter and radiation will be spread ever more thinly and will exert less gravitational influence, but the dark energy will not be diluted; its negative pressure also means that it maintains a constant (or nearly-so) energy density for all time. Therefore, the universe will continue to expand at an accelerating rate, and in another ten billion years it will be almost twice as large as it is now.

I mentioned earlier that the precise value of the $w$ parameter would have a big impact on the universe in the long term. If we run the simulation again with a few different representative values for $w$, this is what we get:

The change in $w$ doesn’t make very much difference on this timescale. You can appreciate why it’s proving difficult for astronomers to get a precise value for $w$ from observation! However, let’s widen our horizons a bit and look at what happens over the next 100 billion years:

Here, even a small change in $w$ makes a dramatic difference. The larger the magnitude of $w$, the faster the expansion accelerates. Sooner or later, we expect the accelerating expansion due to dark energy to tear apart all structures in the universe, starting with the largest (galactic superclusters) and proceeding eventually to the smallest (solid bodies like stars and planets, and eventually even atoms). However, the nature of the dark energy—the as-yet-unknown value of $w$—will determine just how long this process will take.

In a follow-up post, I explore how the size of the cosmic horizon—the region of space that we can observe, given the limitations of lightspeed and the expansion of the universe—changes over time. That will let us get an idea of how long various structures, such as galaxies, might survive.

For reference, and in case you want to mess around with it yourself, here is the Python code I used to generate the data for the preceding graphs.