---

License

DarkEnergyIntro

Wed May 27 16:23:16 2020\ Copyright 2020\ Sandro Dias Pinto Vitenti vitenti@uel.br

---

---

DarkEnergyIntro\ Copyright (C) 2020 Sandro Dias Pinto Vitenti vitenti@uel.br

numcosmo is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

numcosmo is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

---

Dark energy tutorial

In this tutorial we are going to learn how dark energy was necessary to explain the Hubble diagram when the firsts Supernovae Ia (SnIa) data was obtained. In this notebook you will learn:

• What are the basic cosmological parameters.
• How the basic cosmological parameters affect the comoving and luminosity distances.
• How to fit these parameters to adjust the luminosity distance using type Ia supernovae data.

Remember to press shift+enter to execute the cells. All cells must be execute in order!

Loading NumCosmo

The first step is to load both NumCosmo and NumCosmoMath libraries. We also load python libraries necessary to make plots and other miscellaneous tools.

Initializing the library

We now need to initialize the library and redirect its output to this notebook.

Cosmological model - $\Lambda$CDM

We will use the $\Lambda$CDM cosmological model that consists in using the Friedmann metric and filling the universe with radiation, baryons, Cold Dark Matter (CDM) and the cosmological constant ($\Lambda$). In simple terms, the Friedmann metric (with flat spatial surfaces) can be understood as a scale factor $a(t)$ that multiply all rulers (e.g., $x$, $y$ and $z$ directions in Cartesian coordinate system) by $a(t)$ at a time $t$.

Friedmann model

Imagine the following scenario, initially (at time $t_0$) a $30\,\mathrm{cm}$ rubber ruler with no resistance to be stretched is floating in space. After some time $\Delta t$ passes, the size of the rubber ruler is now $30\,\mathrm{cm}\;a(t_0+\Delta t) / a(t_0)$. For instance, if $a(t_0+\Delta t) / a(t_0) = 2$, then at time $t=t_0 + \Delta t$ the ruler is $60\,\mathrm{cm}$ long. However, note that the markings on the ruler still go from $0\,\mathrm{cm} - 30\,\mathrm{cm}$. The absolute value of the difference between two of these markings only matches the physical distance between them at $t_0$. Therefore, in order to measure a real physical distance at any other time $t$, we first measure the coordinates using the ruler and then multiply the result by $a(t) / a(t_0)$. We call these markings as comoving coordinates. In short, the real physical distance $d$ will be always the comoving distance $d_c$ times $a(t)/a(t_0)$, i.e., $d = d_c a(t)/a(t_0)$, and the comoving distances for free falling objects never changes.

Now, let's say we have two non-moving floating balls also in space. We use our rubber ruler and measure their distance as $d(t_0) = d_c(t_0) = 10\,\mathrm{cm}$ at $t_0$. Then, after 2 seconds passes, we measure the distance again at $t = t_0 + 2\,\mathrm{s}$, the rubber ruler indicates again $10\,\mathrm{cm}$ in its markings, so, we need to multiply by $a(t)/a(t_0)$ to get the real distance. Let's say $a(t)/a(t_0) = 4$, consequently, at $t$ the physical distance is $d(t) = 10\,\mathrm{cm}\; a(t) / a (t_0) = 40\,\mathrm{cm}$ while the comoving distance is still $d_c(t) = 10\,\mathrm{cm}$. Since the comoving distance is a constant, we will drop the function argument and label it as $d_c$ only. This computation means that after the time lapse $\Delta t = t - t_0 = 2\,\mathrm{s}$ the distance increased by $30\,\mathrm{cm}$. Here, we should stress that the balls are not moving! What is happening is that the space is stretching and the physical distance is increasing for that reason only. In this example the distance increased by $30\,\mathrm{cm}$ in $2\,\mathrm{s}$, that shows an increase rate of

$$\bar{v}(t) = \frac{d(t_0 + \Delta t) - d(t_0)}{\Delta\,t} = 15\,\mathrm{cm/s}.$$

Dividing this rate by the initial distance $d(t_0) = d_c = 10\,\mathrm{cm}$, we define

$$ \bar{H}(t) \equiv \frac{\bar{v}(t)}{d_c(t)} = \frac{d_c a(t)/a(t_0) - d_c}{\Delta t\,d_c} = \frac{d_c a(t) - d_c a(t_0)}{\Delta t\,d_c a(t_0)} = \frac{1}{a(t_0)}\frac{a(t_0 + \Delta t) - a(t_0)}{\Delta t} = 1.5\,\mathrm{s}^{-1}.$$

The rate above shows that the distance between the balls increase on average by a multiplicative factor $1.5$ each second. Nevertheless, note that there is nothing special about these two balls. Any two non-moving objects will have their initial distance increasing on average by the same factor of $1.5$ per second. This is evinced by the fact that $d_c$ appears in both numerator and denominator in the expression for $\bar{H}(t)$ and, consequently, it cancels out!

The way we defined $\bar{H}(t)$ depends on the elapsed time between the two measurement events $t_0$ and $t = t_0 + \Delta t$. If we consider smaller and smaller time lapses $\Delta t$, we get a better and better approximation for the instantaneous rate the distance increases with time. For this reason we define the instantaneous expansion rate as

$$H(t) \equiv \lim_{\Delta t\to 0}\bar{H}(t).$$

This function is known as Hubble-Lemaître function. When calculated today $t_0$, it is called as Hubble constant $H_0$. The current measurements of $H_0$ indicates a value around $H_0 = 68\,\mathrm{Km/s/Mpc}$, where $\mathrm{Mpc}$ stands for mega parsec, a unit of distance. Let's now call NumCosmo to find out how many meters there are in a $\mathrm{Mpc}$.

The result of the cell above shows that one $\mathrm{pc}$ corresponds to approximately $3.1 \times 10^{16}$ meters or $\approx 3.3$ light-years! To get to $\mathrm{Mpc}$ we still need to multiply by one million, that is, $1 \mathrm{Mpc} \approx 3.1 \times 10^{22}\,\mathrm{m}$. As a reference we also print the value of one parsec in terms of light-years, thus, one parsec is roughly 3 light-years.

Now we can rewrite the Hubble constant in more familiar units.

This is a really small quantity. It means that if you have two floating balls 1 meter apart, they will move away with an initial velocity of $2.2\times 10^{-18}\,\mathrm{m/s}$! Moreover, if there is any kind of resistance, such as, air viscosity, static friction, electromagnetic binding force, etc, the expansion will not provide enough energy to start a movement.

Friedmann equation

The Friedmann model consists in relating the Hubble-Lemaître function $H(t)$, discussed above, with the energy content of the universe, which is encoded in the total energy density $\rho$. The relation useful in our context will be the Friedmann equation

\begin{equation} H^2(t) = \frac{\kappa}{3}\rho(t). \end{equation}

That is, the expansion rate squared is proportional to the energy density of the energy content of the universe. The proportionality factor is

$$\frac{\kappa}{3} = \frac{8\pi\mathrm{G}}{3c^2}.$$

Here $\mathrm{G}$ stands for the Newton gravitational constant and $\mathrm{c}$ is the speed of light. Let's compute it:

The proportionality factor here have units of meter per kilogram. Recall that kinetic energy is mass times velocity squared ($K = mv^2/2$). So, we can rewrite the units above as $\mathrm{m/kg = (m^3/(J s^2))}$, where we used Joule's definition, i.e., $\mathrm{J = kg\, m^2 / s^2}$. We are getting close to something more reasonable now. It seems that the proportionality factor is just one over density but we have an additional factor of one over second squared. However, that's exactly the units on the left hand side of the Friedmann equation, so, dividing both sides by $\kappa/3$, we get

$$\rho_{\mathrm{crit}}(t) \equiv \frac{3H^2(t)}{\kappa} = \rho(t).$$

We are now equating two densities, the left hand side (called critical density) is a simple combination of physical constants and the Hubble function, and the right hand side is the actual energy density of the universe. We can now compute the left hand side today:

That seems to be a very small density. We can compute a more familiar quantity: the mass density. We just divide both sides by $\mathrm{c}^2$. Let's print this density and, for comparison, the density respective to five protons in a cubic meter.

That's roughly the same value, that is, the energy density of the universe computed from the Hubble constant is approximately five protons per square meter, a really small density.

That may look odd, but here we need to remark what we are doing. This density seems small for our standards, nonetheless, this model represents the universe as a whole. So this density must be calculated considering all energy content of the universe and then divide by its size. Here we arrived again in another strange concept, the size of the universe. How could us compute that? To explain this, we need to introduce the concept that is behind this modeling:

The homogeneity and isotropy of the universe.

To build the Friedmann model, we assumed that the universe is homogeneous. This was already evident in our initial description using the floating balls. In that mental exercise, we implicitly assumed that it didn't matter where the balls were, the distance between them would increase no matter their absolute position. Not only that the distance would increase at the same rate, no matter if they were aligned in one direction or another. Strictly speaking, we assumed that the expansion was homogeneous (equal expansion rate independent of the absolute position) and also isotropic (equal expansion rate independent of the direction of the objects considered).

How can this concept be applied in the real world? We can see that our nearby environment is not homogeneous. Clearly the air density is different from the density of a chair. To compute a density we choose a volume, add up all energy content inside of it and then divide this content by the volume. If we choose a volume of a size $\mathrm{cm}^3$ then place it inside a wall, we get the wall density. Then, moving this volume to a empty space, we get the air density. Thus, moving this volume around, we have the density as a function of the position of the volume and will be a rather inhomogeneous space (e.g., a room). We can smooth this function increasing the size of this volume. If the volume is large enough, it will always have some walls, some furniture, air, etc, and the density will not vary so much, but it still will not be homogeneous. Here come the hypothesis, we assume that, if we increase this volume in order to encompass several clusters of galaxies, then at this scale the universe will be homogeneous. That is, it does not matter where the volume is centered, the density will always be the same. Although this is an hypothesis, we can actually measure it! If we survey the night sky and count the number of stars, galaxies, etc, we can compute the mean density and check if it varies as we move the observation point in the sky. We call the scale where the energy distribution gets constant as homogeneity scale.

The energy content

Within the modeling above, the energy content are given in terms of their density. However, now we know that a density needs a distance scale to be defined. For this reason, all densities are assumed to be given in the homogeneity scale such that they are all constants at any time $t$.

Let's start by examining the ordinary matter. Most of it is composed by Baryons (protons, neutrons, etc), which then forms the stars, planets, galaxies, etc. Of course, these objects are in movement (planets orbit stars, stars orbit the galaxy center, etc). Nevertheless, if we consider only the center of mass of the objects inside a volume, it will move less and less as we go to higher scales. This happens because we assume that these movements will cancel out. For this reason, we treat the ordinary matter at the homogeneity scale as free falling dust. But remember, here each grain of dust actually represents a large enough volume such that the region in question is homogeneous! Again, that hypothesis can be tested, we can observe if there is an overall movement of objects even in the largest scales observed, and up to now no such movement was found.

In practice our dust-like fluid behaves as free falling objects. Remember that, as discussed above, the free falling objects do not move with respect to one another. If we have $N$ objects with total mass $M$ inside a volume $V$ at $t_0$, then, at time $t$, we will have the same $N$ objects with total mass $M$ but now in a volume $V[a(t)/a(t_0)]^3$. Thus, the initial density $\rho_b(t_0) = M/V$ becomes

$$\rho_b(t) = \frac{M}{V}\left[\frac{a(t_0)}{a(t)}\right]^3 = \rho_b(t_0)\left[\frac{a(t_0)}{a(t)}\right]^3.$$

In principle, we could apply the same reasoning to the dark matter (DM), however, since we do not see it directly, we do not know if they are in movement. Fortunately, there are other ways to determine if DM is hot, warm or cold (that are outside the scope of this introduction), and observations favor of a cold dark matter CDM. Thus, the density associated to it will behave exactly as the Baryons above, that is,

$$\rho_c(t) = \frac{M_\mathrm{CDM}}{V}\left[\frac{a(t_0)}{a(t)}\right]^3 = \rho_c(t_0)\left[\frac{a(t_0)}{a(t)}\right]^3.$$

Finally, we need to describe the radiation within the universe. It can be composed of light (photons) and also of other ultra-relativistic particles (such as neutrinos). Naturally, they will be moving at the speed of light. However, this movement must be such that the mean amount of radiation inside a given volume will always be the same. Following the reasoning above it may look like it will also behaves like a dust-fluid, and that's almost true. The difference is that the energy associated with radiation depends on the wave-length. But since the universe is expanding, the wave-length is also being stretched. Thus, even though the expanding volume will contain on average the same amount of radiation, however, now the same radiation will have larger wave-length. Since the energy is inversely proportional to the wave-length and all radiation is stretched in the same way, the total energy in a expanding volume $E$ at $t_0$ will be $E[a(t_0)/a(t)]$ at $t$. Then, the radiation energy density will behave as

$$ \rho_r(t) = \frac{E}{V}\left[\frac{a(t_0)}{a(t)}\right]^4 = \rho_r(t_0)\left[\frac{a(t_0)}{a(t)}\right]^4.$$

Note that the radiation energy density drops with the fourth power the scale factor whereas matter drops with the third power.

The cosmological constant $\Lambda$

As its name states, $\Lambda$ density behaves as a constant in both space and time! It is a really strange concept: if we have some density $\rho_\Lambda$ in an initial volume $V$, then in another time $t$ we will have the same amount $\rho_\Lambda$. This means that the total energy associated with this component actually grows when the space expands! Strictly speaking the energy associated with the cosmological constant in a expanding volume initially $V$ at $t_0$ evolves as

$$E_\Lambda(t) = \rho_\Lambda V \left[\frac{a(t)}{a(t_0)}\right]^3.$$

That means that the space expansion actually creates energy.

The dimensionless densities

Now that we have defined all necessary energy content of our model, we need to put them in a convenient form. By manipulating the Friedmann equation, we saw that it can be just the comparison of the actual energy content and the critical density. If we divide both sides of the equation

$$\rho_{\mathrm{crit}}(t) \equiv \frac{3H^2(t)}{\kappa} = \rho(t),$$

by $\rho_{\mathrm{crit}}(t)$, we get

$$\Omega(t) \equiv \frac{\rho(t)}{\rho_{\mathrm{crit}}(t)} = 1.$$

Moreover, writing the total density in terms of the individual ones obtained above results in

$$\frac{\rho_b(t)+\rho_c(t)+\rho_r(t)+\rho_\Lambda}{\rho_{\mathrm{crit}}(t)} = \Omega_b(t) + \Omega_c(t) + \Omega_r(t) + \Omega_\Lambda(t) = 1, \qquad \Omega_i(t) \equiv \frac{\rho_i(t)}{\rho_{\mathrm{crit}}(t)}, \; \mathrm{for\,i\,=\,}b,\,c,\,r,\,\Lambda.$$

We know how each individual density evolves, therefore, we just need to know their value at a given time $t_0$ in order to fix these quantities. Here, we fix their values today, that is, we have four parameters $\Omega_{b0} \equiv \Omega_b(t_0)$, $\Omega_{c0} \equiv \Omega_c(t_0)$, $\Omega_{r0} \equiv \Omega_r(t_0)$ and $\Omega_{\Lambda0} \equiv \Omega_\Lambda(t_0)$. Note that, these parameters are not completely free. Friedmann equation constrain them by imposing that their sum must be equal to one. Thus, in practice, we can just fix three parameters and compute, say, $\Omega_\Lambda$ using the Friedmann equation, that is,

$$\Omega_{\Lambda0} = 1 - \left(\Omega_{b0} + \Omega_{c0} + \Omega_{r0}\right).$$

Creating the cosmological model

We are finally ready to create our cosmological model:

The first line above creates a cosmological model of the type NcHICosmoDEXcdm. That means, a Homogeneous and Isotropic cosmology (HICosmo) with Dark Energy (DE) and cold dark matter (cdm). The X stands for a constant equation of state for the DE, but that's out of the scope of this notebook, so here we just need to fix the $w$ parameter to $w = -1$ to get the behavior of $\Lambda$ described above. The next line changes the parameter from $\Omega_\Lambda to \Omega_k$. This new variable $\Omega_k$ is actually another component of the Friedmann equation that we will not discuss here. In this notebook we just set it to zero.

Let's set the parameters for an universe with no cosmological constant and no radiation. For technical reasons, the radiation density is set using the radiation temperature today $T_{\gamma0}$.

Below we create a model set (MSet), in this simple example it will contain just the HICosmo model, and then we use one of its functionalities to print out the current model parameters.

The model set NcmMSet method above "prettylog" shows all models it has inside and their parameter values. Note that it didn't display $\Omega\Lambda$. The reason is that it is no longer a parameter of the model, since it is completely defined by the Friedmann equation. If we need to see its value we can call the HICosmoDE method "E2Omega_de". This method computes the density as a function of the redshift $z$ and multiplied by

$$E(t) \equiv \frac{H(t)}{H_0},$$

namely $E^2(z)\Omega_\Lambda(z)$.

The call above shows that the amount of cosmological constant is exactly zero. That's because we set $\Omega_{c0} = 0.95$, $\Omega_{b0} = 0.5$ and $\Omega_{r0} = 0$, hence, naturally

$$\Omega_{\Lambda0} = 1 - \left(\Omega_{b0} + \Omega_{c0} + \Omega_{r0}\right) = 0.$$

To plot the densities as they evolve in time, we first define a better "time". As discussed above, the radiation waves inside an expanding universe have their wave-lengths stretched as the universe expand. In the visible light, as we move to larger and larger wave-lengths we get redder and redder light. Thus, we call redshift $z$ the factor defined by

$$z \equiv \frac{a(t_0)}{a(t)} - 1.$$

It is zero when $t = t_0$ and increases as $a(t) < a(t_0)$. Since the universe is expanding, that means that $t < t_0$, that is, $t$ is a moment in the past. That may seem odd but it actually appropriated for cosmology. We are always observing light (and now gravitational waves!) coming from past events. In short, $z = 0$ refers to the present time and larger values to moments in the past.

The code below creates an interactive plot showing the densities as a function of the redshift $z$. You can use the sliders to change the values of the cosmological parameters and see the effect on the densities.

Note that initially $\Omega_\Lambda$ corresponds to roughly 70% of the energy density. Nevertheless, for higher and higher redshifts this trend changes, the density of matter and radiation grow and their contribution quickly surpass the cosmological constant. Remember that, higher redshifts correspond to further in the past and we are in an expanding universe, consequently, as we move to the past the volumes shrink and the densities increase! The radiation density grows only larger than matter's in the far past. Why does that happen?

Cosmological distances

Many important cosmological observables depend on the cosmological distances. So, what are exactly cosmological distances? As we stated before, most of the time we measure light that comes from the sky, that's our main source of information about the universe! Then, it is natural to ask, how long the light that we see now traveled? If we knew the time interval $\Delta t$ between the emission of the light $(t_e)$ and the moment we observed it (now $t_0$), it would be easy to compute the distance if the universe was not expanding, i.e., $c\Delta t$. In an expanding universe, during the time the light travels $c\Delta t$ the universe expands by $a(t_0)/a(t_e)$. Thus, the real physical distance between the light and the point where it originates from is $d \approx a(t_0)c\Delta t / a(t_e)$. However, this is only approximately true, the function $a(t)$ can vary continuously with time and, consequently, the expression is a good approximation when $\Delta t$ is very small. Then, to get a good approximation of the total distance traveled, we need to divide the total time interval in small time slices and sum all contribution. Strictly speaking, we need to integrate over the time interval to obtain the exact result, that is,

$$d(t_0,\,t_e) = a(t_0) \int_{t_e}^{t_0} \frac{c\,\mathrm{d}t'}{a(t')}.$$

We can work a little on this expression to get closer to what we already know about the cosmological model. First, remember that we want to use the redshift $z$ as the time variable. Thus, we can compute our small time interval $\Delta t$ in terms of small $\Delta z$. We can write,

$$\Delta z = z(t+\Delta t) - z(t) = \frac{a(t_0)}{a(t+\Delta t)} - \frac{a(t_0)}{a(t)}.$$

Considering the mean expanding rate $\bar{v}$ as we did before, we can write $a(t+\Delta t) = a(t) + \bar{v}\Delta t$. Substituting this expression back on the equation above results in

$$\Delta z = \frac{a(t_0)}{a(t) + \bar{v}\Delta t} - \frac{a(t_0)}{a(t)} = -a(t_0)\frac{\bar{v}\Delta t}{a(t)\left[a(t) + \bar{v}\Delta t\right]} \approx -a(t_0)\frac{\bar{v}\Delta t}{a^2(t)},$$

where in the last equality we are considering $\Delta t$ so small that $a(t) + \bar{v}\Delta t \approx a(t)$. The first surprise in the expression above is the minus sign, but there is nothing wrong with it. It just means what we already knew: when we move to the past ($\Delta t < 0$) the redshift increases $\Delta z > 0$ and vice-versa. Using this expression and our integral above we get to

$$d(z) \equiv d(z,\,0) = -a(t_0) \int_z^0 \frac{c\,\mathrm{d}z'}{H(z')} = a(t_0) \int_0^z \frac{c\,\mathrm{d}z'}{H(z')} = a(t_0)\frac{c}{H_0} \int_0^z \frac{\mathrm{d}z'}{E(z')}.$$

In the expression above, we used that at time $t_0$ the redshift is zero $z(t_0) = 0$ and that we can use the minus sign to invert the integration order (that only means that we are now adding up $\mathrm{d}z \propto -\Delta z = z(t)-z(t+\Delta t)$ instead of $\mathrm{d}z \propto \Delta t = z(t+\Delta t) - z(t)$). The quantity outside the integral is called Hubble radius $R_H$ and we can compute it using:

Then, the distance traveled by a beam of light emitted at $z$ and observed today ($z=0$) will be $R_H$ times the result of the integral above. The best way to get familiar to these expressions is by plotting its results. But first we need another object "NcDistance". This object is responsible to compute all integrals necessary to evaluate different cosmological distances.

The parameter 5.0 used above just indicates to the object to optimize the distances computations up to $z=5$. But note that it can compute distances up to any redshift.

Playing with the cosmological parameters the user can see that by adding more CDM the same redshift $z$ correspond to shorter distances. It is easy to understand the reason for this. First remember that the expansion rate $H(t)$ is related to the energy content through the Friedmann equation

$$H^2(t) = H_0^2\left[\Omega_{\Lambda0} + \Omega_{b0}(1+z)^3 + \Omega_{c0}(1+z)^3 + \Omega_{r0}(1+z)^4\right],$$

where we wrote the densities in terms of the critical density today $\rho_\mathrm{crit0}$ and the scale factors in terms of the redshift $a(t_0)/a(t) = 1 + z$. We can see in the figure legend that the amount of radiation is very small, even if we increase the temperature, thus, it will not strongly affect the expression above. If we increase the amount of matter, then it means that at larger redshifts the expansion rate was larger (remember $H^2 \propto \Omega_{c0}(1+z)^3$). Consequently, the more we add matter faster the universe expanded in the past. That means that the universe went from $a(t)\to a(t_0)$ in less time and, consequently, the amount of distance traveled is smaller since it is proportional to $c\Delta t$.

Luminosity distance

Suppose that we are observing a star in the sky. What is the relation between its luminosity, when the light is emitted at the star's surface, and the luminosity that we see here? To understand this relation let's first divide the light emitted in spherical waves (at each instant a spherical shell of light leaves the star and spread-out in all directions). Then, we assume that we know the star emits light with luminosity $L_0\;[\mathrm{W}]$ measured in unit of power (Watts). Thus, in a small time interval, a very thin spherical shell of light, with radius $r_0$ and thickness $l_0$, leaves the star with total energy $E_0$ and density $\rho_0$, i.e.,

$$E_0 = \frac{l_0}{c}L_0, \qquad \rho_0 = \frac{E_0}{4\pi l_0 r_0^2} = \frac{L_0}{4\pi c r_0^2}.$$

We already know that the light looses energy when the universe expands. For this reason the energy as the wave arrives here is $E = E_0 / (1 + z)$. Moreover, the thickness of the shell was also stretched, being $l = l_0(1+z)$, and the radius is now $r$ (the distance between us and the star), consequently, the observed energy density is

$$ \rho = \frac{E}{4\pi l r^2 } = \frac{E_0}{4\pi l_0 r^2 (1+z)^2} = \frac{L_0}{4\pi c r^2(1+z)^2}.$$

Measuring the luminosity of the star as it arrives here (for example, through its energy density), we define two quantities,

$$ D_L \equiv \sqrt{\frac{L_0}{4\pi c \rho}}, \qquad d_L(z) \equiv (1+z)d(z),$$

where the first quantity $D_L$ can be directly computed by measuring $\rho$ here on Earth and using our knowledge of $L_0$. The second quantity $d_L(z)$, called luminosity distance, can be compute using our cosmological model, the physical distance $d(z)$ that light travelled between its emission at time $t$ and observation at $t_0$, which depends on the scale factor ratio $a(t_0)/a(t)$. That is, we need to know the redshift $z = a(t_0)/a(t) - 1$!

Measuring redshift

We saw above that to compute the luminosity distance using our cosmological model, we need to know the amount of redshift the light suffered when it traveled from the star to us. The redshift corresponds to how much the wave-length of the light is increased due to the expansion of the space-time. So, if somehow we knew the initial wave-length of the light $\lambda$, then by measuring the wave-length we detect ($\lambda_0$) we could infer the redshift using

$$z = \frac{\lambda_0}{\lambda} - 1.$$

How can we know the wave-length of the light emitted by the star? There are many different physical process going on in a star and each radiates light with different wave-lengths. Luckily, one of these processes is the transition of electrons from different bound states in an atom. These processes are interesting because they can only happen in a discrete manner (they are described by quantum mechanics!). In other words, an electron in a Hydrogen atom can move from one discrete bound state to another, emitting electromagnetic waves (actually photons) with a specific wave-length (that is directly associated with the energy difference between the two bound states).

When we observe the light coming from a star and split it in each different wave-length, we are measuring what we call its spectrum (how much light in each frequency). If there is enough of a given element (say Hydrogen) in a starm we can identify the spectral lines associated to it. These lines are always at the same position in the spectrum - remember, they result from a discrete transition in an atom and the energy difference is always the same. This means that we actually know the $\lambda$ of some spectral lines of a star. Therefore, if we measure its spectrum we can compute $z$!

Type Ia Supernovae

Studying carefully the last sections, one can deduce that, measuring the spectrum of a star (or any other astrophysical object) with enough spectral lines allow us to compute (estimate) its redshift. Given a redshift and a cosmological model, one can compute a cosmological distance, e.g., the luminosity distance. Notwithstanding we do not know what is the Universe's cosmological model. On the contrary, the primary goal in Cosmology is to find out the best cosmological model that describes the Universe, and not to compute distances. Even if we were absolutely sure that the Friedmann model was the correct one, which would be the correct values for its parameters?

The plots above show clearly that varying these parameters result in rather different curves for the cosmological distances. How can we actually deduce the parameters using astrophysical observations? It was showed in the last sections that, if we also knew the luminosity of a star, we could also compute the luminosity distance $D_L$. Furthermore, this distance should match the luminosity distance $d_L(z)$ computed using the model, that is, $d_L(z) = D_L$. This means that, if we can measure the spectrum of a star and infer its luminosity $L_0$, we can vary the cosmological parameters until we get an agreement between $d_L(z)$ and $D_L$.

The astrophysical object that allowed this kind of analysis was not a star, but the explosion of one, the event that we call Type Ia supernova. It was an phenomenological model built for these events that allowed the measurement (inference) of both $z$ and $L$ of a single astrophysical observable.

Model fitting

It seems that we have everything necessary to fix the values of our model parameters. We just need to measure a few supernovae explosions (spectrum and luminosity), compute $z$ and $D_L$ and equate $d_L(z) = D_L$. Right?

No. In practice the problem is more complicated than that. All these measurements can be made, but not perfectly! That is, the $z$ and $D_L$ measurements provide estimates and also their respective errors. For example, if we measure $D_L = x$, how confident we are that $D_L = x$ and not $D_L = x + \delta$ or $D_L = x - \delta$? Of course, this depends on how large $\delta$ is. The error is exactly the measure of how large $\delta$ can be given the measurement imprecisions. Strictly speaking, we associate a probability distribution for $\delta$, for instance, we say $\delta$ is zero on average and has $\sim 32\%$ probability of being $\vert\delta\vert > \sigma$.

This probability 32% is not a random number, it comes from a Gaussian distribution (also known as normal distribution). Assuming $\delta$'s probability distribution is normal, there is $\sim 32\%$ probability of finding it larger than one standard deviation $\sigma$. Moreover, a zero mean normal distribution is completely characterized by $\sigma$, that is, the probability density distribution for $\delta$ is

$$P(\delta) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{\delta^2}{2\sigma^2}}.$$

But, why are we interested in the Gaussian distribution? Why should we assume that the probability for $\delta$ is Gaussian? In practice, we do not know the probability distribution for $\delta$, however, there is an important mathematical theorem that can aid us, the central limit theorem. This theorem asserts that combining different independent random factors result in a single Gaussian distribution for this combination as the number of factors grow. For this reason, assuming that the final error in our measurements results from the combination of different factors (different instrument errors, atmospheric conditions on the observation nights, etc), we can hope that the final distribution is Gaussian (or very well approximated by a Gaussian). There are many techniques to check if the Gaussian distribution is a good approximation for a given data set, or even to try to deduce what is the correct distribution. In this notebook we follow the simplest route, we just assume that the errors are indeed Gaussian.

The discussion above showed that we cannot simple equate $d_L(z)$ to $D_L$ and find the right parameters, we actually have a statistical distribution for $\delta = D_L - d_L(z)$. Let's use now real data to see how this analysis can be made in practice.

The cell below load a set of Type Ia Supernovae as they were analyzed in the earlies 2000. We have a set of redshifts $z_1,\,z_2,\,\dots z_N$ and the [distance modulus] $\mu_{1},\,\mu_{2},\,\dots \mu_{N}$ for $N$ supernovae. Note that we have the distance modulus and not the luminosity distance. Nevertheless, there is a very simple relation between them, its usage is mostly for astronomical reasons,

$$\mu \equiv 5\log_{10}\left(\frac{D_L}{1\mathrm{Mpc}}\right) + 25.$$

The line above shows the data related to the first supernova Ia in the catalog (Nc.DataSNIAId.SIMPLE_LEGACY). Here we introduce a further complication. The empirical model for the type Ia supernovae states that, once their data is fitted their absolute magnitude (which is directly related to their luminosity) are all the same $M_0$, but we actually do not know the value of $M_0$. However, here we will not deal with this point and will just use the sample as it is (it already includes an artificial value of $M_0$).

In this example, it seems that our model is actually a bad explanation for the data. The model and data differ by few $\sigma$. Let's try to improve that and fit the model, we plot below the curve $\mu(z)$ using the model and the observed $\mu_i = 5\log_{10}(D_L / 1\mathrm{Mpc}) + 25$ at $z_i$ and their error bars $\sigma_i$. To guide us we also included in the plot the quantities

$$\bar\delta \equiv \sum_{i=1}^N\left[\frac{\mu(z_i) - \mu_i}{\sigma_i}\right], \qquad \mathrm{Var}(\delta) \equiv \sum_{i=1}^N\left[\frac{\mu(z_i) - \mu_i}{\sigma_i}\right]^2.$$

The first gives the mean value of $\delta$ and the second the sum of the squared differences weighted by the standard deviation. Both quantities have a direct relation with the (likelihood)[https://en.wikipedia.org/wiki/Likelihood_function], that is, the probability of observing the data given the model. One very used technique to infer the parameters that best explain the data is to find the set of parameters that maximize this probability. In this setting this is equivalent to minimize the function $\mathrm{Var}(\delta)$. With that in mind, you can try to modify the parameters below to find the point that minimizes $\mathrm{Var}(\delta)$, this point is known as best fit. Is it possible to find a good fit without dark energy (try to increase $\Omega_{c0}$ until $\Omega_\Lambda$ reaches zero)? If it is not possible, what is the smallest $\mathrm{Var}(\delta)$ we can get without dark energy? Repeat this exercise by varying $H_0$.