Spherical Harmonics - learn to roll your own

It doesn't take much reading about the cosmic microwave background before you butt up against the concept of spherical harmonics (aka those funny looking $Y_l^m(\theta,\phi)$ guys) and if you're like me they threw you for a loop at first glance. Typically introduced in the context of studying anisotropies of the CMB they are essentially an alternate way of representing functions defined on the surface of a sphere. In this case the sphere we're interested in is the sky around the Earth and the function to be represented is the intensity of CMB radiation finally arriving at Earth after billions of years flying through space between photon decoupling during recombination and today. I won't go into too much detail about the physics of the CMB in this post but suffice it to say that these anisotropies are the definitive evidence for the Big Bang and cosmological inflation.

Spherical harmonics can be used to represent any arbitrary function as a (potentially infinite) list of coefficients that weight how much each harmonic contributes. If you're coming from an engineering background this may make more sense when thought of like an extension to the Fourier series where an arbitrary function can be decomposed into a (potentially infinite) weighted sum of sinusoids, each representing a specific frequency component of the original signal. Spherical harmonics are much the same but rather than using one-dimensional sinusoids defined along a circle the decomposition uses two-dimensional functions defined along the surface of a sphere. What is important to grok is that any signal can be totally defined solely in terms of the values of the weights for each harmonic.

To get a better understanding of how exactly this works in practice I wrote a small program to generate image maps to visually see how each unweighted harmonic varies across the surface of the sphere. The first step in doing so was getting a basic handle on the math involved. We'll start with the formal definition of our signal of interest $T(\theta,\phi)$ where $Y_l^m$ is a spherical harmonic (specific to degree $l$ and order $m$) and $a_l^m$ is the weighting factor for that specific harmonic.

$$T(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{+l} a_l^m Y_l^m(\theta,\phi)$$

We're interested in displaying each harmonic individually so we'll need the definition of the spherical harmonic $Y_l^m$ with normalization factor $n_l^m$.

$$Y_l^m(\theta,\phi) = n_l^m P_l^m(cos(\theta)) e^{i m \phi}$$
$$n_l^m = \sqrt{\frac{2 l + 1}{4 \pi} \frac{(l - m)!}{(l + m)!}}$$

Keep in mind $a_l^m$ and $e^{i m \phi}$ are both complex quantities even if the end result we're displaying will just be the real value.

$a_l^m Y_l^m(\theta,\phi) = a_l^m n_l^m P_l^m(cos(\theta)) e^{i m \phi}$
$a_l^m Y_l^m(\theta,\phi) = n_l^m P_l^m(cos(\theta)) a_l^m (cos(m \phi) + i\ sin(m \phi))$

$a_l^m Y_l^m(\theta,\phi) = n_l^m P_l^m(cos(\theta)) (Re(a_l^m) + i\ Imag(a_l^m)) (cos(m \phi) + i\ sin(m \phi))$

$Re(a_l^m Y_l^m(\theta,\phi)) = n_l^m P_l^m(cos(\theta)) (Re(a_l^m) cos(m \phi) - Imag(a_l^m) sin(m \phi))$

For the purposes of the image generation I'm simply using $Re(a_l^m) = Imag(a_l^m) = 1$.

$P_l^m$ is the associated Legendre polynomial. For this the easiest way to deal with the associated Legendre polynomials was to simply precompute them as I wasn't able to find a free C# library for computing them dynamically. Since I'm planning on generating 6000x3000 pixel images that are rectilinear projections of the surface of the sphere I'll need 3000 equally spaced samples. This is relatively simple to generate with numerical computing software such as GNU Octave or MATLAB. The script I used in Octave to generate the values for the first 50 associated Legendre polynomials $P_l^m$ of $cos(\theta)$ and export them to CSVs is below.

for i = 0:49 

csvwrite(strcat(strcat("legendres-",num2str(i)),".csv"), legendre(i, cos(0:pi/3000:pi*2999/3000)));

endfor

Note that this only gets you results for $m >= 0$ but we can extend this with the simple relation below.

$$P_l^{-m} = (-1)^m \frac{(l - m)!}{(l + m)!} P_l^m$$

That's all the math we're using! I wrote up a quick script that reads in the CSV output from the Octave script and uses it to calculate $Re(a_l^m Y_l^m(\theta,\phi))$ for each $l$ where $l \leq 50$ and each $m$ where $-l \leq m \leq l$ linearly sampling $\theta$ between $0$ and $\pi$ 3000 times and $\phi$ between $0$ and $2 \pi$ 6000 times. It then takes the resulting data and outputs an image where the minimum value is blue and the maximum value is red. The entire project is available here: https://github.com/btouellette/SphericalHarmonics/

The whole purpose of doing this was to get a better handle on the concept of spherical harmonics by visual inspection so let's take a look at what we end up with!

$l = 0$, $m = 0$

$l = 1$, $m = 0$

$l = 1$, $m = 1$

$l = 2$, $m = 0$

$l = 2$, $m = 1$

$l = 2$, $m = 2$

$l = 3$, $m = 0$

$l = 3$, $m = 1$

$l = 3$, $m = 2$

$l = 3$, $m = 3$

A quick look at the results shows us that $l$ is related to the number of repetitions on the surface of the sphere and $m$ is related to the "orientation" of the repetitions. I've left out the images for negative $m$ as they aren't particularly illuminating but a quick look will confirm that they are simply the same as the positive $m$ results but with a rotation of $\frac{90°}{m}$ applied.

One last look at these, this time from the perspective of an observer inside the sphere looking outward as we do at the cosmos!

The Friedmann Equation from the beginning

The single most important equation of modern cosmology is the Friedmann equation which relates the expansion of space to the energy density and curvature of the universe.
$$\left(\frac{\dot{a}}{a}\right)^2 = \frac{8 \pi G}{3} \rho - \frac{k c^2}{a^2}$$
It was discovered in 1922 as a consequence of general relativity and accurately describes the Hubble expansion despite being introduced seven years before Hubble discovered his law. It applies in a simplified universe which consists only of a perfect fluid and is isotropic (looks the same no matter which way you look) and homogenous (is the same everywhere). The conditions of isotropy and homogeneity are commonly known as the cosmological principle. While it obviously does not hold on small scales, on large enough distance scales the cosmological principle appears to be true. Averaging over large enough distances things seem to look the same in every direction and at every place. The precise scale at which things begin to look statistically homogenous and the validity of the cosmological principle itself have been debated as larger and larger structures are found but there is still strong evidence in support of it in the incredible homogeneity of the cosmic microwave background.

To speak mathematically about the expansion of the universe we first need to define a scale factor $a(t)$ which is a relation between the size of the universe at time $t$ and now (usually referred to as $t_0$) such that at the current time $a=1$. At a point where the universe is half the size it is today $a = 0.5$. The evolution of the scale factor with time tells us exactly how the universe expands or contracts over time. As a consequence of homogeneity the scale factor can only be a function of time, not position.

We can (and will) derive the Friedmann equation from the Einstein equation but we can arrive at the same result with a much simpler argument. First let us define a comoving coordinate, one which moves with the expansion of the universe such that the coordinate of a point does not change with time. The physical distance between two points in the universe is time-dependent due to the expansion of the universe but the (comoving) coordinate distance between two points remains fixed. By reusing the scale factor to define comoving coordinate $R = \frac{r}{a(t)}$ we can arrive at some simple relations to $R$ for distance, velocity, acceleration, and volume respectively.

$r = a R$
$\dot{r} = \dot{a} R$
$\ddot{r} = \ddot{a} R$
$V = \frac{4}{3} \pi a^3 R^3$

Imagine an expanding spherical shell of mass $m$ at a constant comoving $R$ in a homogenous expanding universe. Using the shell theorem (which generalizes to general relativity) the only gravitational potential energy for the shell comes from the mass interior to the shell. Because the shell is defined by a comoving coordinate the total mass interior to the shell remains constant.
$$M_i = V(t) \cdot \rho(t) = \frac{4}{3} \pi r^3(t) \cdot \rho(t) = \frac{4}{3} \pi a^3(t) R^3 \cdot \rho(t)$$
Conservation of mechanical energy of this spherical shell of mass implies the sum of potential and kinetic energy for the shell remains constant:
$$K(t) + U(t) = E$$
Let us introduce a new parameters $k$ (dimensions of $\frac{1}{length^2}$) such that the energy takes the form $- \frac{1}{2} m k c^2 R^2$. Remember that $R$ for a given expanding mass shell is constant with dimensions $length$.
$$\frac{1}{2} m \dot{r}^2(t) - G \frac{M_i m}{r(t)} = - \frac{1}{2} m k c^2 R^2$$
If the total energy of the system is 0 the mass shell comes to rest at infinity ($k=0$, flat). If the total energy of the shell is negative the shell eventually halts and contracts back in on itself ($k>0$, closed, bounded). If the total energy is positive the shell expands forever ($k<0$, open, unbounded). Replacing $\dot{r}(t)$ and $M_i$ and simplifying:
$$\frac{1}{2} m \dot{a}^2(t) R^2 -  \frac{4}{3} \pi G a^2(t) R^2 \rho(t) m = - \frac{1}{2} m k c^2 R^2$$
$$  \dot{a}^2(t) R^2 -  \frac{8}{3} \pi G a^2(t) R^2 \rho(t) = - k c^2 R^2$$
$$  \left( \frac{\dot{a}(t)}{a(t)} \right)^2 -  \frac{8}{3} \pi G \rho(t) = - \frac{k c^2}{a^2(t)}$$
The final equation correctly does not depend on $r$ or $R$ as it must satisfy our initial assumptions of homogeneity. It must only be a function of parameters which are the same everywhere in the universe and for every potential shell. Any possible changes done to this system that don't violate homogeneity or isotropy won't change the energy of the mass shell so $k$ is constant. As an aside in some situations you will see $k$ interpreted as a dimensionless parameter where $k \in \{+1,0,-1\}$ corresponding to our three cases for the total energy of the shell. Doing so simply requires a reinterpretation of the scale factor $a(t)$ to incorporate some of the aspects of $k$ associated with the curvature of the space and to have units of $length$. I will touch on the curvature aspect later but will maintain the convention of $a(t)$ being dimensionless and $k$ being on a continuum.

While this imaginary universe is certainly evocative this derivation is not exactly precise. What exactly is meant by the density $\rho(t)$? The shell theorem still holds in general relativity but it isn't clear whether this density is simply a matter density as would be appropriate in the Newtonian model or if it encompasses non-Newtonian sources of gravity such as radiation or dark energy. Let's take a step back from this imaginary model and examine things from the perspective of general relativity to attempt to clarify this issue.

The cosmological principle severely restricts the form the metric of the universe may take since the metric must be the same at every point in space and in every direction for a given time. More formally a spacetime is only isotropic and homogeneous if it can be divided up into a family of spacelike three-surfaces which are themselves isotropic and homogeneous. The time coordinate serves to identify individual spacelike three-surfaces with constant time across them. The most generic possible metric of any spacetime can be written as:
$$ds^2 = -\alpha c^2 dt^2 + \sum_{i=1}^{3} \beta_i dx^i dt + dS^2$$
Here $dS^2$ is the line element of space (containing no time components) and $\alpha, \beta_1, \beta_2, \beta_3$ are all scaling constants. The first simplification we can perform is to change to a time coordinate in which $\alpha = 1$ via $dt' = \sqrt{\alpha} dt$. This will also change the values of $\beta_i$ but as they are still just constants we arrive at a metric (where $dt$ and $\beta_i$ are potentially different from the previous metric) with the form:
$$ds^2 = -c^2 dt^2 + \sum_{i=1}^{3} \beta_i dx^i dt + dS^2$$
We know that the four-velocity for points in comoving coordinates must be in a timelike direction only as otherwise isotropy would be violated. Spatial velocity for these points would cause things to look different in the direction of the three-velocity compared to other directions (think Doppler effect). More formally the four-velocity of each point in comoving coordinates is orthogonal to any vector which is tangent to the surface of homogeneity (the spacelike three-surface of constant time). Mathematically this tells us that the four-velocity $u^\alpha = (c, 0, 0, 0)$, a tangent vector is of form $x^\alpha = (0, x^1, x^2, x^3)$, and for all tangent vectors:
$$\textbf{u} \cdot \textbf{x} = g_{\mu \nu} u^\mu x^\nu = g_{t i} u^t x^i = g_{t i} x^i = 0 $$
Since $x^i$ is not necessarily 0 this implies $g_{t i} = 0$ for all $i$. As $g_{i t} = 0$ as well this further implies a lack of mixed space and time terms in the metric so $\beta_i = 0$ resulting in the metric:
$$ds^2 = -c^2 dt^2 + dS^2$$

The three-dimensional metric $dS^2$ is also restricted by assumptions of isotropy and homogeneity. Due to homogeneity no matter what form the metric takes the curvature of space must be the same at all spatial positions for a given time. This limits us to one of three possible configurations. It is either flat, positively curved (spherical), or negatively curved (hyperbolic).

Remember that these images are two-dimensional surfaces embedded in three-dimensions and our use of these terms refers to the three-dimensional universe we live in and not to a lower dimensional embedding. Do not read too heavily into parallels with lower dimensional examples. The important takeaway is that in two-dimensional space these figures are the only types of surfaces which have the same curvature at each point along the two dimensional plane. Similarly if the cosmological principle holds the universe will have the same curvature at each point across the entire three dimensional space. One of the interesting parallels that does hold is that in a positively curved space the sum of the angles of a triangle is greater than 180° while in a negatively curved space the sum will be less than 180°. This is a consequence of the fact that the angle between two lines does not remain constant in a curved space and gives us an avenue for testing the curvature of space. While we can't construct cosmologically sized triangles to measure the angles of we can model the observational effects of different curvature on things like the incoming light from the cosmic microwave background or counts of galaxies at different distances compared to models of galactic formation.

$dS^2$ may vary with time so long as it varies everywhere identically. One source of possible time variation is the scale factor. It must apply everywhere across the space (no position dependence) the same to maintain homogeneity giving $dS^2 = a^2(t) (d\Omega_3)^2$. The other similar source of time dependence is in the curvature which we have not assumed to be constant.

We know that $(d\Omega_3)^2$ for the flat, isotropic, and homogenous universe takes the form $dx^2 + dy^2 + dz^2$ in Cartesian coordinates as it is just the standard metric of three dimensional space. To explore the effect that the curvature might have on the metric we will introduce a fourth dimension $w$ and consider the three dimensional universe as an embedded surface in a four dimensional space. For the spherical case and at a specific timeslice the surface with radius of curvature $R_c$ (not to be confused with our comoving coordinate) is defined by:
$$R_c^2 = x^2 + y^2 + z^2 + w^2$$
We know that the three dimensional radius is defined by:
$$r^2 = x^2 + y^2 + z^2$$
$$w^2 = R_c^2 - r^2$$
Taking the differential of the original equation for the surface:
$$2 R_c\, dR_c = 2 x\, dx + 2 y\, dy + 2 z\, dz + 2 w\, dw$$
Since this slice of geometry is taken at a specific time the curvature is just a constant:
$$0 = x\, dx + y\, dy + z\, dz + w\, dw$$
$$dw = - \frac{x\, dx + y\, dy + z\, dz}{w}$$
$$dw = - \frac{r\, dr}{\sqrt{R_c^2 - r^2}}$$ 

$$dw = - \frac{r\, dr}{\sqrt{R_c^2 - r^2}}$$

We now arrive at the distance function across the space without reference to our higher dimension $w$:

$$(d\Omega_3)^2 = dx^2 + dy^2 + dz^2 + dw^2 = dx^2 + dy^2 + dz^2 + \frac{r^2 dr^2}{R_c^2 - r^2}$$

We recognize $dx^2 + dy^2 + dz^2$ as our standard metric of three dimensional space and we have a an additional term which represents the change in the distance due to the curvature of the space. Rewriting entirely in spherical coordinates results in:

$$(d\Omega_3)^2 = dr^2 + r^2 (d\theta^2 + sin^2\theta\, d\phi^2) + \frac{r^2 dr^2}{R_c^2 - r^2}$$

Combining the $dr^2$ terms and simplifying results in:

$$(d\Omega_3)^2 = \left(1+\frac{r^2}{R_c^2 - r^2}\right) dr^2 + r^2 (d\theta^2 + sin^2\theta\, d\phi^2)$$

$$(d\Omega_3)^2 = \frac{R_c^2}{R_c^2 - r^2} dr^2 + r^2 (d\theta^2 + sin^2\theta\, d\phi^2)$$

$$(d\Omega_3)^2 = \frac{dr^2}{1 - \frac{r^2}{R_c^2}} + r^2 (d\theta^2 + sin^2\theta\, d\phi^2)$$

We can generalize this to the flat case by considering flat as the spherical case but with infinite radius of curvature and we can see as expected the final curvature related term in the metric drops out. For the hyperbolic case we repeat this procedure with a different equation defining the embedding:

$$-R_c^2 = x^2 - y^2 - z^2 - w^2$$
$$-r^2 = x^2 - y^2 - z^2$$
$$w^2 = R_c^2 - r^2$$
$$0 = x\, dx - y\, dy - z\, dz - w\, dw$$

$$dw = \frac{x\, dx - y\, dy - z\, dz}{w}$$
$$dw = -\frac{r\, dr}{\sqrt{R_c^2 - r^2}}$$ 

$$(d\Omega_3)^2 = dx^2 - dy^2 - dz^2 - dw^2 = dx^2 + dy^2 + dz^2 + \frac{r^2 dr^2}{R_c^2 - r^2}$$
$$(d\Omega_3)^2 = dr^2 + r^2 (d\theta^2 + sin^2\theta\, d\phi^2) + \frac{r^2 dr^2}{R_c^2 - r^2}$$

Here we reach the same result as previously except with $R_c^2$ interpreted as the negative of the square of the radius of curvature. We gain some insight into the role of $k$ now by associating it with the Gaussian curvature of a two-dimensional spatial slice of the space at the current time ($a(t) = 1$, necessary since expansion will reduce the curvature) allowing us to express all three cases with a single metric.

$$k = \frac{1}{R_c^2}$$

$$(d\Omega_3)^2 = \frac{dr^2}{1 - k\, r^2} + r^2 (d\theta^2 + sin^2\theta\, d\phi^2)$$

This makes the metric which is the most general possible metric compatible with the cosmological principle:
$$ds^2 = -c^2 dt^2 + a^2(t) \left(\frac{dr^2}{1 - k\, r^2} + r^2 (d\theta^2+sin^2\theta\, d\phi^2)\right)$$
The terms closed/bounded and open/unbounded used for the different cases $k>0$ and $k<0$ refer to the volume integral $V = \int{A\, d\Omega_3}$ where $A$ is constructed taking into account the scale factor. If $k>0$ this integral is definite and the universe can be said to have a finite volume at a fixed time.

Now that we have the final metric to play with we can get really down and dirty with it. Starting with the Einstein equation we look only at the temporal equation:
$$R_{\mu \nu} - \frac{1}{2}g_{\mu \nu}R - \Lambda g_{\mu \nu} = -\frac{8 \pi G}{c^4} T_{\mu \nu}$$
$$R_{t t} - \frac{1}{2}g_{t t}R - \Lambda g_{t t} = -\frac{8 \pi G}{c^4} T_{t t}$$
A quick note here that the stress-energy tensor $T_{\mu \nu}$ will have opposite sign depending on the definition used for the Ricci tensor so this will sometimes appear slightly differently. The sign of the cosmological constant $\Lambda$ varies with the metric signature as well. By the way the history of the cosmological constant is quite interesting. It was originally inserted into the Einstein equation as an effort to describe a static universe and Friedmann and Hubble's discoveries led to it being removed. Recent observations calculating the acceleration of the expansion of the universe have shown the need for an energy density of empty space (or vacuum energy density) which takes the form of the cosmological constant. Moving forward we know $g_{t t} = -c^2$ from a quick look at the metric reducing this to:
$$R_{t t} + \frac{1}{2}c^2 R + \Lambda c^2 = -\frac{8 \pi G}{c^4} T_{t t}$$
The first thing we need to find is the Ricci curvature tensor $R_{\mu \nu}$ which we find by contracting the Riemann tensor as defined by the Christoffel symbols of our metric. We'll use this to get the value of the Ricci scalar $R$ as well:
$$R_{\mu \nu} = R^\alpha_{\mu \alpha \nu} = \partial_\nu \Gamma^\alpha_{\mu \alpha} - \partial_\alpha \Gamma^\alpha_{\mu \nu} + \Gamma^\lambda_{\mu \alpha} \Gamma^\alpha_{\lambda \nu} - \Gamma^\lambda_{\mu \nu} \Gamma^\alpha_{\lambda \alpha}$$
$$R = g^{\mu \nu} R_{\mu \nu}$$
We'll need a bunch of Christoffel symbols so let's go ahead and calculate all of them. This is made much easier by the fact that our metric is diagonal. The only non-zero metric components are matched terms which are listed below along with all their partial derivatives:

$g_{t t} = -c^2$
$g_{r r} = \frac{a^2}{1 - k r^2}$
$g_{\theta \theta} = a^2 r^2$
$g_{\phi \phi} = a^2 r^2 sin^2\theta$

$\partial_t g_{t t} = 0$

$\partial_r g_{t t} = 0$

$\partial_\theta g_{t t} = 0$

$\partial_\phi g_{t t} = 0$

$\partial_t g_{r r} = \frac{2 a \dot{a}}{1 - k r^2}$

$\partial_r g_{r r} = \frac{2 k r a^2}{(1 - k r^2)^2}$

$\partial_\theta g_{r r} = 0$

$\partial_\phi g_{r r} = 0$

$\partial_t g_{\theta \theta} = 2 a \dot{a} r^2$

$\partial_r g_{\theta \theta} = 2 a^2 r$

$\partial_\theta g_{\theta \theta} = 0$

$\partial_\phi g_{\theta \theta} = 0$

$\partial_t g_{\phi \phi} = 2 a \dot{a} r^2 sin^2\theta$

$\partial_r g_{\phi \phi} = 2 a^2 r sin^2\theta$

$\partial_\theta g_{\phi \phi} = 2 a^2 r^2 sin\theta cos\theta$

$\partial_\phi g_{\phi \phi} = 0$

Since our metric is diagonal this also makes finding the contravariant components trivial using $g_{a a}g^{a a} = 1$:

$g^{t t} = -\frac{1}{c^2}$
$g^{r r} = \frac{1 - k r^2}{a^2}$
$g^{\theta \theta} = \frac{1}{a^2 r^2}$
$g^{\phi \phi} = \frac{1}{a^2 r^2 sin^2\theta}$

Now we have everything we need to calculate the Christoffel symbols using:

$$\Gamma^i_{k \ell} = \frac{1}{2} g^{i m}( \partial_\ell g_{m k} + \partial_k g_{m \ell} - \partial_m g_{k \ell})$$

The metric connection is necessarily torsion-free $\Gamma^i_{k \ell} = \Gamma^i_{\ell k}$ reducing the number of calculations necessary.

$\Gamma^t_{t t} = \frac{1}{2} g^{t m}( \partial_t g_{m t} + \partial_t g_{m t} - \partial_m g_{t t}) = 0$

$\Gamma^t_{t r} = \Gamma^t_{r t} = \frac{1}{2} g^{t m}( \partial_r g_{m t} + \partial_t g_{m r} - \partial_m g_{t r}) = 0$

$\Gamma^t_{t \theta} = \Gamma^t_{\theta t} = \frac{1}{2} g^{t m}( \partial_\theta g_{m t} + \partial_t g_{m \theta} - \partial_m g_{t \theta}) = 0$

$\Gamma^t_{t \phi} = \Gamma^t_{\phi t} = \frac{1}{2} g^{t m}( \partial_\phi g_{m t} + \partial_t g_{m \phi} - \partial_m g_{t \phi}) = 0$

$\Gamma^t_{r r} =  \frac{1}{2} g^{t m}( \partial_r g_{m r} + \partial_r g_{m r} - \partial_m g_{r r}) = -\frac{1}{2} g^{t t} \partial_t g_{r r} = \frac{a \dot{a}}{c^2 (1 - k r^2)}$

$\Gamma^t_{r \theta} = \Gamma^t_{\theta r} = \frac{1}{2} g^{t m}( \partial_\theta g_{m r} + \partial_r g_{m \theta} - \partial_m g_{r \theta}) = 0$

$\Gamma^t_{r \phi} = \Gamma^t_{\phi r} = \frac{1}{2} g^{t m}( \partial_\phi g_{m r} + \partial_r g_{m \phi} - \partial_m g_{r \phi}) = 0$

$\Gamma^t_{\theta \theta} = \frac{1}{2} g^{t m}( \partial_\theta g_{m \theta} + \partial_\theta g_{m \theta} - \partial_m g_{\theta \theta}) = - \frac{1}{2} g^{t t} \partial_t g_{\theta \theta} = \frac{a \dot{a} r^2}{c^2}$
$\Gamma^t_{\theta \phi} = \Gamma^t_{\phi \theta} = \frac{1}{2} g^{t m}( \partial_\phi g_{m \theta} + \partial_\theta g_{m \phi} - \partial_m g_{\theta \phi}) = 0$

$\Gamma^t_{\phi \phi} = \frac{1}{2} g^{t m}( \partial_\phi g_{m \phi} + \partial_\phi g_{m \phi} - \partial_m g_{\phi \phi}) = - \frac{1}{2} g^{t t} \partial_t g_{\phi \phi} = \frac{a \dot{a} r^2 sin^2\theta}{c^2}$

$\Gamma^r_{t t} = \frac{1}{2} g^{r m}( \partial_t g_{m t} + \partial_t g_{m t} - \partial_m g_{t t}) = 0$

$\Gamma^r_{t r} = \Gamma^r_{r t} = \frac{1}{2} g^{r m}( \partial_r g_{m t} + \partial_t g_{m r} - \partial_m g_{t r}) = \frac{1}{2} g^{r r} \partial_t g_{r r} = \frac{\dot{a}}{a}$

$\Gamma^r_{t \theta} = \Gamma^r_{\theta t} = \frac{1}{2} g^{r m}( \partial_\theta g_{m t} + \partial_t g_{m \theta} - \partial_m g_{t \theta}) = 0$

$\Gamma^r_{t \phi} = \Gamma^r_{\phi t} = \frac{1}{2} g^{r m}( \partial_\phi g_{m t} + \partial_t g_{m \phi} - \partial_m g_{t \phi}) = 0$

$\Gamma^r_{r r} =  \frac{1}{2} g^{r m}( \partial_r g_{m r} + \partial_r g_{m r} - \partial_m g_{r r}) = \frac{1}{2} g^{r r} \partial_r g_{r r} = \frac{k r}{(1 - k r^2)}$

$\Gamma^r_{r \theta} = \Gamma^r_{\theta r} = \frac{1}{2} g^{r m}( \partial_\theta g_{m r} + \partial_r g_{m \theta} - \partial_m g_{r \theta}) = 0$

$\Gamma^r_{r \phi} = \Gamma^r_{\phi r} = \frac{1}{2} g^{r m}( \partial_\phi g_{m r} + \partial_r g_{m \phi} - \partial_m g_{r \phi}) = 0$

$\Gamma^r_{\theta \theta} = \frac{1}{2} g^{r m}( \partial_\theta g_{m \theta} + \partial_\theta g_{m \theta} - \partial_m g_{\theta \theta}) = - \frac{1}{2} g^{r r} \partial_r g_{\theta \theta} = - r (1 - k r^2)$
$\Gamma^r_{\theta \phi} = \Gamma^r_{\phi \theta} = \frac{1}{2} g^{r m}( \partial_\phi g_{m \theta} + \partial_\theta g_{m \phi} - \partial_m g_{\theta \phi}) = 0$

$\Gamma^r_{\phi \phi} = \frac{1}{2} g^{r m}( \partial_\phi g_{m \phi} + \partial_\phi g_{m \phi} - \partial_m g_{\phi \phi}) = - \frac{1}{2} g^{r r} \partial_r g_{\phi \phi} = - r sin^2\theta (1 - k r^2)$

$\Gamma^\theta_{t t} = \frac{1}{2} g^{\theta m}( \partial_t g_{m t} + \partial_t g_{m t} - \partial_m g_{t t}) = 0$

$\Gamma^\theta_{t r} = \Gamma^\theta_{r t} = \frac{1}{2} g^{\theta m}( \partial_r g_{m t} + \partial_t g_{m r} - \partial_m g_{t r}) = 0$

$\Gamma^\theta_{t \theta} = \Gamma^\theta_{\theta t} = \frac{1}{2} g^{\theta m}( \partial_\theta g_{m t} + \partial_t g_{m \theta} - \partial_m g_{t \theta}) = \frac{1}{2} g^{\theta \theta} \partial_t g_{\theta \theta} = \frac{\dot{a}}{a}$

$\Gamma^\theta_{t \phi} = \Gamma^\theta_{\phi t} = \frac{1}{2} g^{\theta m}( \partial_\phi g_{m t} + \partial_t g_{m \phi} - \partial_m g_{t \phi}) = 0$

$\Gamma^\theta_{r r} =  \frac{1}{2} g^{\theta m}( \partial_r g_{m r} + \partial_r g_{m r} - \partial_m g_{r r}) = 0$

$\Gamma^\theta_{r \theta} = \Gamma^\theta_{\theta r} = \frac{1}{2} g^{\theta m}( \partial_\theta g_{m r} + \partial_r g_{m \theta} - \partial_m g_{r \theta}) = \frac{1}{2} g^{\theta \theta} \partial_r g_{\theta \theta} = \frac{1}{r}$

$\Gamma^\theta_{r \phi} = \Gamma^\theta_{\phi r} = \frac{1}{2} g^{\theta m}( \partial_\phi g_{m r} + \partial_r g_{m \phi} - \partial_m g_{r \phi}) = 0$

$\Gamma^\theta_{\theta \theta} = \frac{1}{2} g^{\theta m}( \partial_\theta g_{m \theta} + \partial_\theta g_{m \theta} - \partial_m g_{\theta \theta}) = 0$
$\Gamma^\theta_{\theta \phi} = \Gamma^\theta_{\phi \theta} = \frac{1}{2} g^{\theta m}( \partial_\phi g_{m \theta} + \partial_\theta g_{m \phi} - \partial_m g_{\theta \phi}) = 0$

$\Gamma^\theta_{\phi \phi} = \frac{1}{2} g^{\theta m}( \partial_\phi g_{m \phi} + \partial_\phi g_{m \phi} - \partial_m g_{\phi \phi}) = - \frac{1}{2} g^{\theta \theta} \partial_\theta g_{\phi \phi} = - sin\theta cos\theta$

$\Gamma^\phi_{t t} = \frac{1}{2} g^{\phi m}( \partial_t g_{m t} + \partial_t g_{m t} - \partial_m g_{t t}) = 0$
$\Gamma^\phi_{t r} = \Gamma^\phi_{r t} = \frac{1}{2} g^{\phi m}( \partial_r g_{m t} + \partial_t g_{m r} - \partial_m g_{t r}) = 0$
$\Gamma^\phi_{t \theta} = \Gamma^\phi_{\theta t} = \frac{1}{2} g^{\phi m}( \partial_\theta g_{m t} + \partial_t g_{m \theta} - \partial_m g_{t \theta}) = 0$
$\Gamma^\phi_{t \phi} = \Gamma^\phi_{\phi t} = \frac{1}{2} g^{\phi m}( \partial_\phi g_{m t} + \partial_t g_{m \phi} - \partial_m g_{t \phi}) = \frac{1}{2} g^{\phi \phi} \partial_t g_{\phi \phi} = \frac{\dot{a}}{a}$
$\Gamma^\phi_{r r} =  \frac{1}{2} g^{\phi m}( \partial_r g_{m r} + \partial_r g_{m r} - \partial_m g_{r r}) = 0$
$\Gamma^\phi_{r \theta} = \Gamma^\phi_{\theta r} = \frac{1}{2} g^{\phi m}( \partial_\theta g_{m r} + \partial_r g_{m \theta} - \partial_m g_{r \theta}) = 0$
$\Gamma^\phi_{r \phi} = \Gamma^\phi_{\phi r} = \frac{1}{2} g^{\phi m}( \partial_\phi g_{m r} + \partial_r g_{m \phi} - \partial_m g_{r \phi}) = \frac{1}{2} g^{\phi \phi} \partial_r g_{\phi \phi} = \frac{1}{r}$
$\Gamma^\phi_{\theta \theta} = \frac{1}{2} g^{\phi m}( \partial_\theta g_{m \theta} + \partial_\theta g_{m \theta} - \partial_m g_{\theta \theta}) = 0$
$\Gamma^\phi_{\theta \phi} = \Gamma^\phi_{\phi \theta} = \frac{1}{2} g^{\phi m}( \partial_\phi g_{m \theta} + \partial_\theta g_{m \phi} - \partial_m g_{\theta \phi}) = \frac{1}{2} g^{\phi \phi} \partial_\theta g_{\phi \phi} =  \frac{cos\theta}{sin\theta}$
$\Gamma^\phi_{\phi \phi} = \frac{1}{2} g^{\phi m}( \partial_\phi g_{m \phi} + \partial_\phi g_{m \phi} - \partial_m g_{\phi \phi}) = 0$

Great, now we have everything we need to find $R_{\mu \nu}$ using the contraction of the Riemann tensor and the following observations which help to reduce the number of necessary calculations

1. The Ricci tensor is symmetric on a Riemannian manifold ($R_{\mu \nu} = R_{\nu \mu}$)
2. If the first and second or second and third lower index in the Riemann tensor are matched that component of the tensor must be zero through the skew symmetry of the tensor
3. The Christoffel symbols of this metric are always zero if they do not have exactly two of their three indices matched
4. If an upper and lower index of the Christoffel symbol are the matching ones it is only non-zero if the non-matching index precedes the matching one in this list: $t, r, \theta, \phi$
5. If the Christoffel symbol has a $t$ lower index it is only a function of $t$, not any spatial coordinate

$$R_{\mu \nu} = R^\alpha_{\mu \alpha \nu} = R^t_{\mu t \nu} + R^r_{\mu r \nu} + R^\theta_{\mu \theta \nu} + R^\phi_{\mu \phi \nu}$$
$R_{t t} = R^t_{t t t} + R^r_{t r t} + R^\theta_{t \theta t} + R^\phi_{t \phi t}$
$R_{t t} = R^r_{t r t} + R^\theta_{t \theta t} + R^\phi_{t \phi t}$

$R_{t t} = \partial_t \Gamma^r_{t r} + \Gamma^r_{t r} \Gamma^r_{r t} + \partial_t \Gamma^\theta_{t \theta} + \Gamma^\theta_{t \theta} \Gamma^\theta_{\theta t} + \partial_t \Gamma^\phi_{t \phi} + \Gamma^\phi_{t \phi} \Gamma^\phi_{\phi t}$
$R_{t t} = 3 \partial_t \left(\frac{\dot{a}}{a}\right) + 3 \frac{\dot{a}^2}{a^2}$
$R_{t t} = 3 \left(\frac{\ddot{a}}{a} - \frac{\dot{a}^2}{a^2} \right) + 3 \frac{\dot{a}^2}{a^2}$
$R_{t t} = 3 \frac{\ddot{a}}{a}$

$R_{t r} = R_{r t} = R^t_{t t r} + R^r_{t r r} + R^\theta_{t \theta r} + R^\phi_{t \phi r}$
$R_{t r} = R_{r t} = R^\theta_{t \theta r} + R^\phi_{t \phi r}$

$R_{t r} = R_{r t} = \Gamma^\theta_{t \theta} \Gamma^\theta_{\theta r} - \Gamma^r_{t r} \Gamma^\theta_{r \theta} + \Gamma^\phi_{t \phi} \Gamma^\phi_{\phi r} - \Gamma^r_{t r} \Gamma^\phi_{r \phi} = 0$

$R_{t \theta} = R_{\theta t} = R^t_{t t \theta} + R^r_{t r \theta} + R^\theta_{t \theta \theta} + R^\phi_{t \phi \theta}$
$R_{t \theta} = R_{\theta t} = R^r_{t r \theta} + R^\phi_{t \phi \theta}$

$R_{t \theta} = R_{\theta t} = \Gamma^r_{t r} \Gamma^r_{r \theta} - \Gamma^\theta_{t \theta} \Gamma^r_{\theta r} + \Gamma^\phi_{t \phi} \Gamma^\phi_{\phi \theta} - \Gamma^\theta_{t \theta} \Gamma^\phi_{\theta \phi} = 0$

$R_{t \phi} = R_{\phi t} = R^t_{t t \phi} + R^r_{t r \phi} + R^\theta_{t \theta \phi} + R^\phi_{t \phi \phi}$
$R_{t \phi} = R_{\phi t} = R^r_{t r \phi} + R^\theta_{t \theta \phi}$

$R_{t \phi} = R_{\phi t} = \Gamma^r_{t r} \Gamma^r_{r \phi} - \Gamma^\phi_{t \phi} \Gamma^r_{\phi r} + \Gamma^\theta_{t \theta} \Gamma^\theta_{\theta \phi} - \Gamma^\phi_{t \phi} \Gamma^\theta_{\phi \theta} = 0$

$R_{r r} = R^t_{r t r} + R^r_{r r r} + R^\theta_{r \theta r} + R^\phi_{r \phi r}$
$R_{r r} = R^t_{r t r} + R^\theta_{r \theta r} + R^\phi_{r \phi r}$

$R_{r r} = - \partial_t \Gamma^t_{r r} + \Gamma^r_{r t} \Gamma^t_{r r} + \partial_r \Gamma^\theta_{r \theta} + \Gamma^\theta_{r \theta} \Gamma^\theta_{\theta r} - \Gamma^t_{r r} \Gamma^\theta_{t \theta} - \Gamma^r_{r r} \Gamma^\theta_{r \theta}+ \partial_r \Gamma^\phi_{r \phi} + \Gamma^\phi_{r \phi} \Gamma^\phi_{\phi r} - \Gamma^t_{r r} \Gamma^\phi_{t \phi} - \Gamma^r_{r r} \Gamma^\phi_{r \phi}$

$R_{r r} = - \partial_t \left( \frac{a \dot{a}}{c^2 (1 - k r^2)} \right) + \frac{\dot{a}}{a} \frac{a \dot{a}}{c^2 (1 - k r^2)} + \partial_r \left( \frac{1}{r} \right) + \frac{1}{r^2} - \frac{a \dot{a}}{c^2 (1 - k r^2)} \frac{\dot{a}}{a} - \frac{k r}{1 - k r^2} \frac{1}{r} + \partial_r \left( \frac{1}{r} \right) + \frac{1}{r^2} - \frac{a \dot{a}}{c^2 (1 - k r^2)} \frac{\dot{a}}{a} - \frac{k r}{1 - k r^2} \frac{1}{r}$

$R_{r r} = - \frac{\dot{a}^2}{c^2 (1 - k r^2)} - \frac{a \ddot{a}}{c^2 (1 - k r^2)} - \frac{2 k r}{1 - k r^2} \frac{1}{r} - \frac{a \dot{a}}{c^2 (1 - k r^2)} \frac{\dot{a}}{a}$

$R_{r r} = - \frac{2 \dot{a}^2}{c^2 (1 - k r^2)} - \frac{a \ddot{a}}{c^2 (1 - k r^2)} - \frac{2 k}{1 - k r^2}$

$R_{r \theta} = R_{\theta r} = R^t_{r t \theta} + R^r_{r r \theta} + R^\theta_{r \theta \theta} + R^\phi_{r \phi \theta}$
$R_{r \theta} = R_{\theta r} = R^t_{r t \theta} + R^\phi_{r \phi \theta}$

$R_{r \theta} = R_{\theta r} = \Gamma^\phi_{r \phi} \Gamma^\phi_{\phi \theta} - \Gamma^\theta_{r \theta} \Gamma^\phi_{\theta \phi} = 0$

$R_{r \phi} = R_{\phi r} = R^t_{r t \phi} + R^r_{r r \phi} + R^\theta_{r \theta \phi} + R^\phi_{r \phi \phi}$
$R_{r \phi} = R_{\phi r} = R^t_{r t \phi} + R^\theta_{r \theta \phi}$

$R_{r \phi} = R_{\phi r} = \Gamma^\theta_{r \theta} \Gamma^\theta_{\theta \phi} - \Gamma^\phi_{r \phi} \Gamma^\theta_{\phi \theta} = 0$

$R_{\theta \theta} = R^t_{\theta t \theta} + R^r_{\theta r \theta} + R^\theta_{\theta \theta \theta} + R^\phi_{\theta \phi \theta}$
$R_{\theta \theta} = R^t_{\theta t \theta} + R^r_{\theta r \theta} + R^\phi_{\theta \phi \theta}$

$R_{\theta \theta} = - \partial_t \Gamma^t_{\theta \theta} + \Gamma^\theta_{\theta t} \Gamma^t_{\theta \theta} - \partial_r \Gamma^r_{\theta \theta} + \Gamma^\theta_{\theta r} \Gamma^r_{\theta \theta} - \Gamma^t_{\theta \theta} \Gamma^r_{t r} - \Gamma^r_{\theta \theta} \Gamma^r_{r r} + \partial_\theta \Gamma^\phi_{\theta \phi} + \Gamma^\phi_{\theta \phi} \Gamma^\phi_{\phi \theta} - \Gamma^t_{\theta \theta} \Gamma^\phi_{t \phi} - \Gamma^r_{\theta \theta} \Gamma^\phi_{r \phi}$

$R_{\theta \theta} = - \partial_t \left( \frac{a \dot{a} r^2}{c^2} \right) + \frac{\dot{a}}{a} \frac{a \dot{a} r^2}{c^2} + \partial_r \left( r (1 - k r^2) \right) - \frac{1}{r} r (1 - k r^2) - \frac{a \dot{a} r^2}{c^2} \frac{\dot{a}}{a} + r (1 - k r^2) \frac{k r}{(1 - k r^2)} + \partial_\theta \left( \frac{cos\theta}{sin\theta} \right) + \frac{cos\theta}{sin\theta} \frac{cos\theta}{sin\theta} - \frac{a \dot{a} r^2}{c^2} \frac{\dot{a}}{a} + r (1 - k r^2) \frac{1}{r}$

$R_{\theta \theta} = - \left( \frac{\dot{a}^2 r^2}{c^2} + \frac{a \ddot{a} r^2}{c^2} \right) + \left( (1 - k r^2) - 2 k r^2 \right) + k r^2 - \left( \frac{1}{sin^2\theta} \right) + \frac{cos^2\theta}{sin^2\theta} - \frac{\dot{a}^2 r^2}{c^2}$

$R_{\theta \theta} = - \frac{a \ddot{a} r^2}{c^2} + (1 - 2 k r^2) - \left( \frac{1}{sin^2\theta} \right) + \frac{cos^2\theta}{sin^2\theta} - \frac{2 \dot{a}^2 r^2}{c^2}$

$R_{\theta \theta} = - \frac{a \ddot{a} r^2}{c^2} + (1 - 2 k r^2) + \frac{cos^2\theta - 1}{sin^2\theta} - \frac{2 \dot{a}^2 r^2}{c^2}$

$R_{\theta \theta} = - \frac{a \ddot{a} r^2}{c^2} - 2 k r^2 - \frac{2 \dot{a}^2 r^2}{c^2}$

$R_{\theta \phi} = R_{\phi \theta} = R^t_{\theta t \phi} + R^r_{\theta r \phi} + R^\theta_{\theta \theta \phi} + R^\phi_{\theta \phi \phi}$
$R_{\theta \phi} = R_{\phi \theta} = R^t_{\theta t \phi} + R^r_{\theta r \phi}$

$R_{\theta \phi} = R_{\phi \theta} = \Gamma^\phi_{\theta \phi} \Gamma^t_{\phi t} - \Gamma^\phi_{\theta \phi} \Gamma^r_{\phi r} = 0$

$R_{\phi \phi} = R^t_{\phi t \phi} + R^r_{\phi r \phi} + R^\theta_{\phi \theta \phi} + R^\phi_{\phi \phi \phi}$
$R_{\phi \phi} = R^t_{\phi t \phi} + R^r_{\phi r \phi} + R^\theta_{\phi \theta \phi}$
$R_{\phi \phi} = - \partial_t \Gamma^t_{\phi \phi} + \Gamma^\phi_{\phi t} \Gamma^t_{\phi \phi} - \partial_r \Gamma^r_{\phi \phi} + \Gamma^\phi_{\phi r} \Gamma^r_{\phi \phi} - \Gamma^t_{\phi \phi} \Gamma^r_{t r} - \Gamma^r_{\phi \phi} \Gamma^r_{r r} - \partial_\theta \Gamma^\theta_{\phi \phi} + \Gamma^\phi_{\phi \theta} \Gamma^\theta_{\phi \phi} - \Gamma^t_{\phi \phi} \Gamma^\theta_{t \theta} - \Gamma^r_{\phi \phi} \Gamma^\theta_{r \theta}$

$R_{\phi \phi} = - \partial_t \left( \frac{a \dot{a} r^2 sin^2\theta}{c^2} \right) + \frac{\dot{a}}{a} \frac{a \dot{a} r^2 sin^2\theta}{c^2} - \partial_r (- r sin^2\theta (1 - k r^2)) + \frac{1}{r} (- r sin^2\theta (1 - k r^2)) - \frac{a \dot{a} r^2 sin^2\theta}{c^2} \frac{\dot{a}}{a} - (- r sin^2\theta (1 - k r^2)) \frac{k r}{(1 - k r^2)} - \partial_\theta (- sin\theta cos\theta) + \frac{cos\theta}{sin\theta} (- sin\theta cos\theta) - \frac{a \dot{a} r^2 sin^2\theta}{c^2} \frac{\dot{a}}{a} - (- r sin^2\theta (1 - k r^2)) \frac{1}{r}$

$R_{\phi \phi} = - \partial_t \left( \frac{a \dot{a} r^2 sin^2\theta}{c^2} \right) + \partial_r (r sin^2\theta (1 - k r^2)) + k r^2 sin^2\theta + \partial_\theta (sin\theta cos\theta) - cos^2\theta - \frac{\dot{a}^2 r^2 sin^2\theta}{c^2}$
$R_{\phi \phi} = - (\frac{\dot{a}^2 r^2 sin^2\theta}{c^2} + \frac{a \ddot{a} r^2 sin^2\theta}{c^2}) + (sin^2\theta (1 - k r^2) - 2 k r^2 sin^2\theta) + k r^2 sin^2\theta + (cos^2\theta - sin^2\theta) - cos^2\theta - \frac{\dot{a}^2 r^2 sin^2\theta}{c^2}$
$R_{\phi \phi} = - \frac{2 \dot{a}^2 r^2 sin^2\theta}{c^2} - \frac{a \ddot{a} r^2 sin^2\theta}{c^2} - 2 k r^2 sin^2\theta $

Whew...

We have the Ricci tensor! It is diagonal just like our metric. We could have saved some time not calculating off diagonal components but they are trivial for the most part and a diagonal metric is not a guarantee that the Ricci tensor is diagonal as well. Luckily getting the Ricci scalar from the Ricci tensor is a lot easier than getting the Ricci tensor from the Riemann tensor.
$$R = g^{\mu \nu} R_{\mu \nu}$$
$$R = -\frac{1}{c^2} 3 \frac{\ddot{a}}{a} - \frac{1 - k r^2}{a^2} \left( \frac{2 \dot{a}^2}{c^2 (1 - k r^2)} + \frac{a \ddot{a}}{c^2 (1 - k r^2)} + \frac{2 k}{1 - k r^2} \right) - \frac{1}{a^2 r^2} \left( \frac{a \ddot{a} r^2}{c^2} + 2 k r^2 + \frac{2 \dot{a}^2 r^2}{c^2} \right) - \frac{1}{a^2 r^2 sin^2\theta} \left( \frac{2 \dot{a}^2 r^2 sin^2\theta}{c^2} + \frac{a \ddot{a} r^2 sin^2\theta}{c^2} + 2 k r^2 sin^2\theta \right)$$
$$R = - \frac{3\ddot{a}}{a c^2} - \left( \frac{2 \dot{a}^2}{a^2 c^2} + \frac{\ddot{a}}{a c^2} + \frac{2 k}{a^2} \right) - \left( \frac{\ddot{a}}{a c^2} + \frac{2 k}{a^2} + \frac{2 \dot{a}^2}{a^2 c^2} \right) - \left( \frac{2 \dot{a}^2}{a^2 c^2} + \frac{\ddot{a}}{a c^2} + \frac{2 k}{a^2} \right)$$
$$R = - \frac{3\ddot{a}}{a c^2} - 3 \left( \frac{2 \dot{a}^2}{a^2 c^2} + \frac{\ddot{a}}{a c^2} + \frac{2 k}{a^2} \right)$$

$$R = - \frac{6\ddot{a}}{a c^2} - \frac{6 \dot{a}^2}{a^2 c^2} - \frac{6 k}{a^2}$$

The only term in the Einstein equation we're missing now is the energy-momentum tensor $T_{\mu \nu}$ (sometimes called the stress-energy tensor). For a perfect fluid in equilibrium (which our cosmological fluid satisfies due to the cosmological principle) this tensor takes a fairly simple form:

$$T_{\mu \nu} = \left( \rho c^2 + p \right) u_\mu u_\nu + p g_{\mu \nu}$$

As we elaborated earlier $u_\alpha = (-c, 0, 0, 0)$ in our comoving coordinates and we can see that the tensor is also diagonal making this:

$T_{t t} = \left( \rho c^2 + p \right) u_t u_t + p g_{t t} =  \left( \rho c^2 + p \right) c^2 - p c^2 = \rho c^4$
$T_{r r} = p g_{r r} = \frac{p a^2}{1 - k r^2}$
$T_{\theta \theta} = p g_{\theta \theta}$
$T_{\phi \phi} = p g_{\phi \phi}$

And we're pretty much done! Swap all the values we've found into the temporal Einstein equation and we get
$$R_{t t} + \frac{1}{2}c^2 R + \Lambda c^2 = -\frac{8 \pi G}{c^4} T_{t t}$$
$$3 \frac{\ddot{a}}{a} + \frac{1}{2}c^2 \left(- \frac{6\ddot{a}}{a c^2} - \frac{6 \dot{a}^2}{a^2 c^2} - \frac{6 k}{a^2} \right) + \Lambda c^2 = -\frac{8 \pi G}{c^4} \rho c^4$$
$$- \frac{3 \dot{a}^2}{a^2} - \frac{3 k c^2}{a^2} + \Lambda c^2 = -8 \pi G \rho$$

$$\left( \frac{\dot{a}}{a} \right)^2 = \frac{8 \pi G}{3} \rho - \frac{k c^2}{a^2} + \frac{\Lambda c^2}{3}$$

This looks (thankfully) quite familiar and we now can identify $\rho$ as the energy density and not just the mass density and we have the addition of the $\Lambda$ term which represents the contribution of a cosmological constant (the energy density of empty space) to the expansion of the universe. In Newtonian mechanics this value is zero however recent research points to the existence of dark energy of which the energy of the vacuum ($\Lambda$) is the simplest possible form.

Itching for more? We've only examined the temporal Einstein equation. There are three more we can calculate for the matched spatial coordinates. All three end up reducing to the second Friedmann equation which describes the second derivative of the scale factor so let's just go with the $r$ equation.

$$R_{r r} - \frac{1}{2} g_{r r} R - \Lambda g_{r r} = -\frac{8 \pi G}{c^4} T_{r r}$$

$$- \frac{2 \dot{a}^2}{c^2 (1 - k r^2)} - \frac{a \ddot{a}}{c^2 (1 - k r^2)} - \frac{2 k}{1 - k r^2} - \frac{1}{2} \frac{a^2}{1 - k r^2} \left( - \frac{6\ddot{a}}{a c^2} - \frac{6 \dot{a}^2}{a^2 c^2} - \frac{6 k}{a^2} \right) - \Lambda \frac{a^2}{1 - k r^2} = -\frac{8 \pi G}{c^4} \frac{p a^2}{1 - k r^2}$$

$$- \frac{2 \dot{a}^2}{c^2} - \frac{a \ddot{a}}{c^2} - 2 k - \frac{a^2}{2} \left( - \frac{6\ddot{a}}{a c^2} - \frac{6 \dot{a}^2}{a^2 c^2} - \frac{6 k}{a^2} \right) - \Lambda a^2 = -\frac{8 \pi G}{c^4} p a^2$$

$$- \frac{2 \dot{a}^2}{a^2 c^2} - \frac{2 k}{a^2} + \frac{2\ddot{a}}{a c^2} + \frac{3 \dot{a}^2}{a^2 c^2} + \frac{3 k}{a^2} - \Lambda = -\frac{8 \pi G}{c^4} p$$

$$\frac{2 \ddot{a}}{a} + \frac{\dot{a}^2}{a^2} + \frac{k c^2}{a^2} - \Lambda c^2 = -\frac{8 \pi G}{c^2} p$$

We're close but we want the resulting equation to be a function of $\frac{\ddot{a}}{a}$ only and not $\frac{\dot{a}}{a}$ so let's substitute the first Friedman equation in.

$$\frac{2 \ddot{a}}{a} + \frac{8 \pi G}{3} \rho - \frac{k c^2}{a^2} + \frac{\Lambda c^2}{3} + \frac{k c^2}{a^2} - \Lambda c^2 = -\frac{8 \pi G}{c^2} p$$

$$\frac{2 \ddot{a}}{a} + \frac{8 \pi G}{3} \rho - \frac{2 \Lambda c^2}{3} = -\frac{8 \pi G}{c^2} p$$

$$\frac{\ddot{a}}{a} = -\frac{4 \pi G}{c^2} p - \frac{4 \pi G}{3} \rho + \frac{\Lambda c^2}{3}$$

$$\frac{\ddot{a}}{a} = - \frac{4 \pi G}{3} \left( \frac{3 p}{c^2} + \rho \right)  + \frac{\Lambda c^2}{3}$$

So we've arrived at the only equations that describe the expansion or contraction of a universe compatible with the cosmological principle. As we can see these are only dependent on the energy density $\rho$, the curvature of space $k$, the cosmological constant aka the vacuum energy $\Lambda$, and the pressure of the cosmological fluid $p$. Different types of energy scale with the size of the universe at different rates and have different pressures associated with them. While the true history of the universe is a constant interplay between the different types we can roughly divide the history of the universe into different eras where these equations are dominated by each type, the radiation, matter, and dark energy eras. Building on top of these equations and bringing in quantum mechanics and astronomical observations has allowed cosmology to develop as a rich and exciting branch of physics.

Configuring MathJax for Blogger

The first step to enabling TeX formatting on your Blogger page is to include the MathJax Javascript file in your Blogger template. Under Template -> Edit HTML you should add the following just after the opening <head> tag:




    <script src="//cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript">    

    MathJax.Hub.Config({

        HTML: ['input/TeX','output/HTML-CSS'],

        TeX: { extensions: ['AMSmath.js','AMSsymbols.js'], equationNumbers: { autoNumber: 'AMS' } },

        extensions: ['tex2jax.js'],

        jax: ['input/TeX','output/HTML-CSS'],

        tex2jax: { inlineMath: [ ['$','$'], ['\\(','\\)'] ],

                   displayMath: [ ['$$','$$'], ['\\[','\\]'] ],

                   skipTags: ["script","noscript","style","textarea","code"],

                   processEscapes: true },

        'HTML-CSS': { availableFonts: ['TeX'], linebreaks: { automatic: true } }

    });

    </script>

In this case I have removed 'pre' from the list of HTML tags skipped over by MathJax. If you are using the Compose tab to post in Blogger rather than raw HTML this is necessary as it uses the <pre> tag throughout. If you want to show TeX intermingled with <code> blocks you can remove that tag from skipTags as well.

MathJax will use the color defined under Page -> Text Color in your template customization so make sure that is set correctly as the Blogger editor will explicitly set the text color in some cases.

This should be all you need to do to enable TeX processing in all posts and comments on your blog. I have listed a few examples along with their raw text below for you to test with.

\(P(E) = {n \choose k} p^k (1-p)^{n-k}\)
\[P(E) = {n \choose k} p^k (1-p)^{n-k}\]

When $a \ne 0$, there are two solutions to \(ax^2 + bx + c = 0\) and they are $$x = {-b \pm \sqrt{b^2-4ac} \over 2a}.$$


\(P(E) = {n \choose k} p^k (1-p)^{n-k}\)
\[P(E) = {n \choose k} p^k (1-p)^{n-k}\]

When $a \ne 0$, there are two solutions to \(ax^2 + bx + c = 0\) and they are $$x = {-b \pm \sqrt{b^2-4ac} \over 2a}.$$