Constraining Cosmological Parameters with Type Ia Supernovae Part I

In this problem, we will use data from Type Ia supernovae to constrain cosmological parameters.

This problem is based on The Supernova Cosmology Project, which is one of the two research teams that first determined a positive cosmological constant using redshift data from Type Ia supernovae. The team took data taken from various telescopes, including the HST to measure the magnitude of various supernovae as a function of their redshift. This work led to the 1997 discovery of the fact that the expansion of the universe is actually accelerating. The same work eventually won the 2011 Nobel Prize and the 2015 Breakthrough Prize in Fundamental Physics.

We use data from the said project for our analysis. The dataset we will use is called the SCP Union2.1 compilation. The compilation has data from 833 supernovae, drawn from 19 datasets but only a certain number pass usability cuts so we are going to use a subset of those 833 supernovae.

You can download the dataset from Download the dataset that says Union2.1 Compilation Magnitude vs. Redshift Table (for your own cosmology fitter). Study the dataset and answer the following question.

  1. How many usable supernovae are in the SCP Union2.1 dataset?

The columns are: supernova name, redshift, distance modulus, distance modulus error, and the probability that the supernova was hosted by a low-mass galaxy. Using this data, we can make what's called a Hubble diagram, which is a plot of distance against redshift.

Make the Hubble diagram with the given dataset. Your plot should look like the following.

The brightness of an object is measured by the apparent magnitude written in terms of the flux, $$m = -2.5 \log_{10}{\left( \frac{f}{f_x} \right)}$$ where $f_x$ is some normalization. In a similar fashion, the absolute magnitude is, $$M = -2.5 \log_{10}{\left( \frac{L}{L_x} \right)}$$ where $L_x$ is also some normalization. $f_x$ and $L_x$ are chosen such that the absolute magnitude corresponds to the apparent magnitude for a body 10 pc away. This makes the equation, $$m - M = 5 \log_{10}{\left( \frac{d_L}{10~\mathrm{pc}} \right)}$$ $$m= M + 5 \log_{10}{\left( \frac{d_L}{1~\mathrm{pc}} \right)} + 25$$ We can then write the magnitude $m$ as, $$ m = 5 \log_{10}\left( H_o d_L \right) + \mathcal{M} $$ $\mathcal{M}$ is a nuisance parameter that encodes uncertainties in the intrinsic luminosity and the Hubble constant. $d_L(z)$ depends on the underlying cosmology and hence, in turn, $m$ depends on the cosmology. This means, if we have a measure of $m$, we can compare it to theoretical prediction of $m$ to test whether or not our universe supports the cosmology. $$m_{\mathrm{th}}(z, \left\{ p_j\right\}) = 5 \log_{10}{\left( \mathrm{H}_o \; d_L(z, \left\{ p_j \right\}) \right)} + \mathcal{M}$$ where $z$ is the redshift and $\left\{ p_j \right\}$ is the set of underlying cosmological parameters. The subscript th just means that this is the expected theoretical magnitude. Recall that $d_L(z)$ depends on the cosmological parameters as, $$ d_L(z) = \frac{H_o^{-1}}{\sqrt{\left| \Omega_K \right|} } \mathrm{sinn} \left[ \sqrt{\left| \Omega_K \right|} \int_0^z \frac{\mathrm{d}z'}{\sqrt{ \Omega_M (1 + z')^3 + \Omega_\Lambda (1 + z')^{3(1 + w)} + \Omega_R (1 + z')^4 + \Omega_K(1 + z')^2 }} \right] $$ where "sinn(x)" is $\sinh{x}$, $\sin{x}$, or $x$ for an open, closed, and flat universe respectively, and $\Omega_K = 1 - \Omega_K - \Omega_R - \Omega_\Lambda$. This means $d_L(z)$ depends on the density parameters $\Omega$ and equation of state $w$. This dependence obviously then holds for $m$, allowing us to use the theoretical vs experimental values of $m$ to constrain those cosmological parameters. We will use maximum likelihood analysis to constrain our parameters. We can use $\chi^2$ value between data and theory to write our likelihood. In particular, $$ \chi^2 = ( {\bf m} - {\bf m}_{\mathrm{th}})^{\mathrm{T}} {\bf \mathrm{Cov}}^{-1} ( {\bf m} - {\bf m}_{\mathrm{th}}) $$ where $\bf \mathrm{Cov}$ is the covariance matrix. Then, the likelihood is just $$ \mathcal{L} = k \exp{\left( - \frac{\chi^2}{2} \right)} $$ where $k$ is a normalization factor. The normalization isn't super important but can be used to ensure that $\mathcal{L}$ is bounded by $1$, which would allow an interpretation in terms of probability. Work out the normalization factor.

  1. What is $k$? Round off to nearest 3 decimal places.
We can now write a Python function for our likelihood. Make sure to handle when your $\chi^2$ value might be nan.

def calculate_likelihood(k, chi2):
    # Implement the function here.
    return likelihood

We're going to need the theoretical values of magnitude $m_{\mathrm{th}}$ so let's implement it but before that we need to implement the luminosity distance $d_L$.

c = 299792  # Speed of light in km/s
H0 = 70 # Hubble constant
M = 25 # Nuisance parameter

def luminosity_distance(z, w, Omega_M, Omega_L):
    # Implement the function here.
    return d_L

def m_th(z, w, Omega_M, Omega_L):
    # Implement the function here. 
    # Calculate d_L with the function luminosity_distance() above.
    return likelihood