In-Class Assignment 25

In-Class Assignment 25#

Due by the end of the day, Wednesday 16 April, 2025

Estimating \(H_0\) from Type Ia Supernovae#

Learning Objectives#

Explore the Phillips Relation using observational data.
Practice plotting data with error bars and fitting to observational data.
Using results of fitted data to compute cosmological parameters.

Type Ia supernova are used as standardizable candles to measure cosmological distances. By observing the lightcurve and measuring how long it takes for the supernova to dim we can empirically determine its brightness via the Phillips relation.

The paper Spectra and Hubble Space Telescope Light Curves of Six Type Ia Supernovae at 0.511 < z < 1.12 and the Union2 Compilation by Amanullah et al. made a data set available that has ~ 500 Type Ia supernovae.

Note

That paper does a far, far more sophisticated analysis than we do, and they fit for other cosmological parameters that we will, so we will not get the same answer as they do. But this is a fun dataset to try out regression with.

import numpy as np
import matplotlib.pyplot as plt

The dataset has 4 columns:

supernova identifier
redshift, \(z\), due to cosmological expansion
distance modulus, \(\mu\), defined as:

\[\mu = m - M = 5 \log_{10} \left (\frac{d}{10~\mbox{pc}}\right )\]

where \(m\) is the apparent magnitude of the supernova (what we observe) and \(M\) is the absolute magnitiude of the supernova (inferred empirically from the lightcurve).
uncertainty in \(\mu\)

a. - Sort and plot the data#

Download the source data for this example here.

First, sort the data by redshift z.
Pass the indexing to create new arrays using the sorted indices.
Plot the distance modulus (\(\mu\)) vs. redshift (\(z\)) and the associated error.

Comment qualitatively on the observed trend in this plot.

#data = np.genfromtxt("SCPUnion2_mu_vs_z.txt",
#                     dtype=[("name", "S6"), ("z", "f8"), ("mu", "f8"), ("dmu", "f8")])

# a result here

Now let’s plot the distance modulus vs redshift

# a plot here

You can compare this result to Figure 9 from the paper mentioned above.

b. - Producing a Hubble diagram#

We can also make a plot that looks a bit more like a Hubble diagram by plotting redshift vs. distance, computing distance as:

\[d = 10 \cdot 10^{\mu/5}~\mathrm{pc}\]

Compute the distance using the above equation.
Plot the redshift as a function of the computed distance (in Gpc) using a scatter plot.

Comment qualitatively on the scatter at large distances.

# b result here

This shows a linear relation between distance and redshift at low redshift.

Cosmological parameters#

We have the distance modulus, which is related to the magnitudes via:

\[\mu = m - M = 5 \log_{10} \left (\frac{d}{10~\mbox{pc}}\right )\]

Since cosmologists usually work in terms of Mpc, let’s rewrite this as:

\[\mu = m - M = 5 \log_{10} \left (\frac{d}{10~\mbox{pc}}\right ) + 5 \log_{10} \left ( \frac{1~\mbox{Mpc}}{1~\mbox{Mpc}}\right ) = 5 \log_{10} \left (\frac{d}{1~\mbox{Mpc}}\right ) + 25\]

Now, in an expanding Universe, the distance that goes here is the luminosity distance which can be expressed via an expansion in redshift as (for \(z \ll 1\)):

\[d_L \approx \frac{c}{H_0} \left [ z + \frac{1}{2} (1 - q_0) z^2 + \ldots \right ]\]

Here \(H_0\) is the Hubble constant and \(q_0\) is the deceleration parameter.

We’ll try to estimate \(H_0\) from this data.

This is different than what the original paper did—they used a value of \(H_0\) to find a \(\Omega_m\) and \(\Omega_\Lambda\) using the more general expression for luminosity distance given in Perlmutter et al. 1997 (see footnote 14):

\[d_L(z; \Omega_M, \Omega_\Lambda, H_0) = \frac{c(1+z)}{H_0 \sqrt{|\kappa|}} S\left ( \sqrt{|\kappa|} \int_0^z \left [(1 + z')^2(1 + \Omega_M z') - z'(2 + z')\Omega_\Lambda\right ]^{-1/2} dz' \right )\]

with \(S(x)\) and \(\kappa\) defined as:

\(\Omega_M + \Omega_\Lambda > 1\): \(S(x) = \sin(x)\); \(\kappa = 1 - \Omega_M - \Omega_\Lambda\)
\(\Omega_M + \Omega_\Lambda = 1\): \(S(x) = x\); \(\kappa = 1\)
\(\Omega_M + \Omega_\Lambda < 1\): \(S(x) = \sinh(x)\); \(\kappa = 1 - \Omega_M - \Omega_\Lambda\)

Using this definition would require us to combine integration and fitting.

c. - Fitting low redshift data to measure \(H_{0}\)#

We want to fit:

\[\mu = 5\log_{10} \left (\frac{cz}{H_0 \cdot 1~\mbox{Mpc}} \left [1 + \frac{1}{2} (1 - q_0) z \right ] \right ) + 25\]

which we’ll write as:

\[\mu = 5\log_{10} \left (a_0 z \left [1 + \frac{1}{2} (1 - a_1) z \right ] \right ) + 25\]

This is a nonlinear expression in terms of the fit parameters, \(a_0\), \(a_1\). Once we get \(a_0\), we can get Hubble’s constant as:

\[H_0 = \frac{c}{a_0 \cdot 1~\mbox{Mpc}}\]

First, use a mask to filter our data to only include data for \(z < 0.2\).
Plot the low redshift data as we did in part a.

# c result here

# c result here

d. - Fitting to the low redshift data#

We’ll use the SciPy fitting routine via method of Least Squares for this. Let’s start by writing the residual that the fit will use.

def resid(avec, z, mu, dmu):
    return (mu - (5 * np.log10(avec[0] * z * (1 + 0.5 * (1 - avec[1]) * z)) + 25)) / dmu

import scipy optimize - from scipy import optimize
define the above residual
Pass the inital guess data to the scipy routine, including the errors on mu.
Use the [0] result from the fit to find \(H_0\).
If there’s time, plot the raw low redshift together, and overlay the estimated non-linear equation from c. and using the values returned from the fit.

Does your computed value agree with the paper?

#from scipy import optimize

#def resid(avec, z, mu, dmu):
#    return (mu - (5 * np.log10(avec[0] * z * (1 + 0.5 * (1 - avec[1]) * z)) + 25)) / dmu

Now let’s start with some initial guesses

# imagine H0 = 50
#c = 3.e5   # km/s 
#H0_guess = 50   # km/s/Mpc

#a0 = c / H0_guess
#a1 = 1

# d result here

From this fit, we can recover the Hubble constant:

# H_0 computed here

Finally, let’s plot our fit

# plot raw data and fit data

We should find a reasonable value of \(H_0 = 69.6~\mathrm{km/s/Mpc}\).