In [ ]:
from __future__ import division, print_function, absolute_import, unicode_literals
import numpy as np
from astropy.table import Table
import matplotlib.pyplot as plt
%matplotlib inline
Today we learned about a variety of different explosive, extragalactic transients. While the lectures focused on recently discovered, rare transients, the most famous transients are without question Type Ia supernovae (SNe). SNe Ia have nearly uniform peak luminosities, which can be standardized via first and second order corrections, such that they are standardizable candels. They are the best distance indicators at high redshift, and the 2011 Nobel Prize was awarded to Adam Riess, Brian Schmidt, and Saul Perlmutter for the discovery of the accelerating universe via Type Ia SNe.
During this exercise we will use PTF data to calibrate the brightness of Type Ia SNe, and then measure Hubble's constant, $H_0$, using this calibration. Following that we will use PTF observations to limit the size of an exploding white dwarf star.
By AA Miller (c) 2016 Jun 14
Previously we learned about image subtraction and the importance of removing the flux from "static" sources (i.e. host galaxies), when trying to measure the brightness of SNe. The PTF public data releases do not include image subtraction products, and the software to perform image subtraction is not (at the moment, anyway) easily installed and implemented in python. Thus, for this exercise we will bypass the public PTF data and instead use published light curves from Firth et al. 2015. Firth et al. (2015) present a study of the rise time of Type Ia SNe discovered by PTF. As they have already performed image subtraction, we will utilize the light curves produced in that study for our exercise. Table 2 of Firth et al. contains the light curves, and that data can be accessed in data/Firth14Tbl2.txt
. One brief note before we start, that study includes SNe from the La Silla-QUEST survey, which we have commented out in the file in the data/
directory. We do this to ensure all light curves were taken in the same filter.
As a first step we will examine the light curves and the formatting of the data file. As before, we will read the data into an astropy
table file.
In [ ]:
# execute this cell
SNlcs = Table.read("../data/Firth14Tbl2.txt", format = 'ascii')
SNlcs
The first thing to notice is that the table does not include magnitude measurements. Gaaaasp The horror!!
As an important point of background - Firth et al. normalize all images to a common zero-point of 27.0 mag (AB), and then after image subtraction perform forced PSF photometry at the location of the PTF SN on the subtraction image. One advantage of this method is that it enables the measurement of negative fluxes (which for SNe isn't particularly useful, but for variable stars is extremely important for detecting events such as eclipses).
For SNe, which should have no flux in the reference image (though there are rare cases where this may not be the case), to calculate magnitudes from flux (or counts
in the case of the Firth et al. study) use the following equation:
where $m$ is the magnitude, $\textrm{ZP}$ is the zero-point, and $f$ is the flux (or counts).
As a first step, convert the counts measurements for each SN to magnitudes.
Problem A1 Convert counts in the SN light curves table to magnitudes using the equation provided above. Store the results in an array called mag
.
In [ ]:
mag = # complete
The previous line should have led to a NumPy run error: and thus, we have encountered one of the downsides of using negative flux measurements - it is not possible to take the $\log$ of a negative number.
Furthermore, the magnitude array that we just created includes measurements where the signal-to-noise ratio (SNR) is $< 1$. Typically, astronomical sources are only considered detected when their flux exceeds some threshold, usually defined in units of the noise (e.g., 3$\sigma$, 5$\sigma$, or in very conservative cases 10$\sigma$).
Problem A2 Create a boolean array, called det
, which tracks epochs in which the SNe are actually detected.
Hint - you must choose the limits at which the SN is considered detected.
In [ ]:
det = # complete
In addition to tracking the mag of the SNe at each epoch, we need to track the uncertainty on each of those measurements. To convert uncertainties in flux to mag uncertainties use the following equation:
$$\Delta m = 1.0857 \frac{\sigma_f}{f},$$where $\Delta m$ is the mag uncertainty and $\sigma_f$ is the flux uncertainty. [You can arrive at this equation by differentiating the first equation in this notebook.]
Problem A3 Calculate the uncertainties for the magnitude measures at each epoch in the SN light curve table and store the results in an array called mag_unc
.
In [ ]:
mag_unc = # complete
Now that we have calculated magnitudes, let's examine some light curves. As discussed earlier today, an important aspect of SN science is determining whether or not a new discovery is a young SN. Examining non-detections, and the corresponding flux limits, can constrain the age of a SN at the time of discovery [more on this later]. Thus, a careful measurement of the upper limits is very important. As an example, plot the light curve for the first SN in the table, PTF09dsy, which is a typical SN from Firth et al.
Problem B1 Plot the light curve, including uncertainties, of PTF09dsy. For epochs where the SN is not detected, plot upper limits (use the v
symbol in matplotlib).
In [ ]:
# complete
plt.errorbar( # complete
plt.ylim( # complete
plt.xlabel( # complete
plt.ylabel( # complete
Problem B2 Did PTF discover PTF09dsy at an early epoch?
Type your response to B2 here
Now examine the family of PTF light curves to see how they compare.
Problem B3 In a single figure, plot all of the PTF light curves from Firth et al. For this figure ignore upper limits. Label each light curve so you can tell them apart via a legend.
Hint - use a for loop to keep your code clean and simple. It may also be useful to increase the size of this particular figure for clarity.
In [ ]:
plt.figure( # complete
# complete
plt.plot( # complete
# complete
plt.legend( # complete
Problem B4 Which SN stands out above the rest from the Firth et al. sample?
Type your response to B4 here
Before we can use Type Ia SNe to measure distances (and, eventually, $H_0$) we must first calibrate their luminosities. Virtually all distance indicators are calibrated in the same way via the "distance ladder." Geometric parallax determines the distance to pulsating variables (e.g., RR Lyrae stars, Cepheids), calibrating those stars, which are then used to calibrate other distance indicators in nearby galaxies. [Note - historically there are many many rungs on the distance ladder which involve many steps for calibration. Recent work, however, has directly calibrated Cepheids via parallax leading to a two step calibration Cepheids $\rightarrow$ SNe Ia, e.g. Riess et al. (2016).]
At the end of problem 1, you should have identified PTF11kly as the especially unique SN in the Firth et al. sample. And, indeed, PTF11kly (aka SN 2011fe) is unique in many, many respects. PTF detected this SN shortly after explosion (more on that later...) in the spectacular Pinwheel Galaxy, M 101. Several hundred papers have been written on this SN that was discovered less than 5 years ago! M 101 is close enough to the the Milky Way that it is possible to detect individual Cepheids with the Hubble Space Telescope, which enables a precise distance measurement via the Cepheid Period-Luminosity relation. Thus, we can use PTF11kly to calibrate the absolute magnitude of Type Ia SNe, thereby measuring distances to all the PTF SNe in the Firth et al. sample.
In addition to being significantly brighter than the other SNe in the sample, the PTF11kly light curve is different for another reason as well. The following is a bit different from the other problems we have encountered so far in that it is a bit open ended, however, it is an important exercise in data exploration.
Problem A1 Identify the difference, aside from peak brightness, between PTF11kly and the other light curves in the Firth et al. sample.
Hint - don't spend too much time on this as the answer is below, but also don't just scroll down without trying to figure this out first.
In [ ]:
# demonstrate that PTF11kly is different from the other SNe
Earlier it was mentioned that Type Ia SNe are standardizable, but there was no discussion of how they are standardizable. For now we will assume that all Type Ia SNe have the same absolute magnitude (a proxy for luminosity for our 1 filter light curves) at peak. In practice, detailed light-curve fitting algorithms such as SALT or MLCS are used to standardize the luminosities of SNe Ia, but the use of these tools is beyond the scope of this problem. Furthermore, SALT, MLCS, and all precise SN distance measurement techniques require light curves in at least two filters, which is not available for PTF data.
To calibrate the absolute magnitude of PTF11kly, we need to determine the distance to M 101, in units of mag, a.k.a the distance modulus.
Problem B1 Determine the distance to M 101, and store the result in mu_M101
.
Hint - you might find the answer on NED or in Riess et al. (2016).
In [ ]:
mu_M101 =
Above, you determined that PTF11kly is different from the other SNe in that PTF observed it in the $g$-Band. SNe do not have flat spectra (in an AB sense) so we cannot calibrate the absolute magnitude in one band and apply the results to another filter. Thus, we will introduce our first cheat of the workshop, which is that we will use a non-PTF light curve to calibrate PTF11kly in the $R$-band. [In addition to being the wrong filter, PTF also did not observe PTF11kly at peak.]
Fortunately, the KAIT telescope observed PTF11kly in the $R$-band covering the peak of the SN light curve. [The KAIT $R$ filter and PTF $R$ filter are not identical, but we will ignore those differences for now. The KAIT light curve shows that PTF11kly peaked at $R = 10.02 \; \mathrm{mag}$.
Problem B2 Store the peak $R$-band brightness in a variable called peak11kly
.
In [ ]:
peak11kly =
Problem B3 Determine the peak absolute magnitude of Type Ia SNe in the $R$-band, store the result in a variable called M_Ia
. Confirm that your result makes sense.
In [ ]:
M_Ia =
Now that we have calibrated the peak absolute magnitude of Type Ia SNe, we can measure $H_0$. Prior to measuring $H_0$, we will try to develop some intuition for the uniformity of SNe Ia at peak.
Traditionally, after SNe candidates are discovered they are sent to the IAU for confirmation, after which they are officially named SN YYYY??, where YYYY is the year, and ?? is an alphabetical sequence following the order in which the SNe were discovered. [Note - with modern surveys discovering hundreds of new SNe every year this scheme is no longer used.] It has been said that one can make a low-scatter Hubble diagram using the SN redshift and mag at discovery from the IAU Circulars, without any sort of filter corrections.
We can test this hypothesis using the PTF data in Firth et al. If SNe are standard candles, and the hypothesis is true, then a plot of $\mathrm{mag}_\mathrm{disc}$ vs. $z$ should show small scatter.
Problem A1 Plot $R$ at the epoch of discovery against redshift for the SNe in the Firth et al. sample. Store the results in arrays called disc_mag
, disc_mag_unc
, and z
.
In [ ]:
# complete
# complete
disc_mag = # complete
disc_mag_unc = # complete
z = # complete
# complete
# complete
plt.errorbar( # complete
plt.ylim(# complete
plt.xlabel(# complete
plt.ylabel(# complete
Based on your above plot, are SNe (at the time of discovery) good standard candles?
While the correlation in the above plot is weak, and the scatter large, we maintain that a decent Hubble diagram can be made using the information present in old IAU circulars.
Problem A2 List some reasons the above exercise would not work well for PTF, but it could work well for previous surveys.
Hint - think back to Mansi's talk from this morning.
Type your response to A2 here
Of course, we know SNe Ia are standard(-izable) candles [they wouldn't award a Nobel Prize for something's that wrong, right? ... Right?!], and that they are ~standard at the time of peak. Now we will see if PTF SNe Ia are standard candles near the time of peak.
Problem A3 Plot the peak $R$ mag for each of the SNe as a function of redshift. Store the results in arrays called mag_peak
and mag_peak_unc
.
Hint - no need for fancy fitting, simply use the brightest observation of each SNe. If you are looking for a challenge you can fit the light curves and interpolate to get the peak, however, note that there is no simple functional form to fit, so you may get worse results following this procedure.
In [ ]:
mag_peak = # complete
mag_peak_unc = # complete
# complete
# complete
plt.errorbar(# complete
plt.ylim(# complete
plt.xlabel(# complete
plt.ylabel(# complete
How does the correlation and scatter look now?
Now that we have empirically demonstrated a correlation between peak brightness and distance for SNe Ia [this is not exactly true, but let's roll with it], we can use the fact that they are standard candles to infer the distance to each. This will allow us to determine distance as a function of recession velocity, aka Hubble's Constant.
Problem B1 Using PTF11kly as a calibrator, determine the distance, in Mpc, to each of the other SNe in the Firth et al. sample. Recall that the distance modulus, $\mu$, is given by:
$$\mu = m - M = 5 \log_{10} (\frac{d}{10}),$$where $m$ is the observed mag, $M$ the absolute mag, and $d$ the distance in pc.
In [ ]:
d = # complete
Problem B2 Plot recession velocity as a function of distance for the PTF Type Ia SNe, thus making a version of the Hubble diagram. How do your results compare to Hubble's original diagram?
In [ ]:
# complete
plt.errorbar(# complete
plt.xlabel(# complete
plt.ylabel(# complete
Problem B3 Perform a linear-least squares fit to the data in the previous plot to determine the value of $H_0$ from PTF SNe. Then, plot the line corresponding to the best fit on the previous plot.
Hint - there are many ways to perform linear-least squares in Python including:
np.polyfit()
np.linalg.lstsq()
scipy.stats.linregress()
and, of course, for a problem this simple it would also be straight-forward to directly code the result yourself.
Note For an actual publication, performing an ordinary least squares (OLS) fit to this data would be inappropriate as the distance measurements have significant uncertainties. Furthermore, flipping the axes, such that the fit is d vs. v, and inverting the slope to get $H_0$ also would not work, as fitting $X$ vs. $Y$ is different from fitting $Y$ vs. $X$. For a brief tutorial on a better approach in this specific case, see Hogg et al. 2010.
In [ ]:
# complete
print('The fit results in H_0 = {:.1f} km/s/Mpc'.format(# complete
Riess et al. (2016) use a large sample of Cepheids and SNe to measure $H_0 = 73.24 \pm 1.74$. [Some of you may be familiar with the tension in $H_0$ measurements between SNe and the Planck measurements of the cosmic microwave background. This is a time-domain workshop, so we are definitely #TeamSNe.]
Problem B4 How does the Riess et al. measurement compare to what you derived in the previous problem?
Type your response to B4 here
Not bad for a basic method and only 10 SNe!!
Now we will significantly change pace, as we pivot away from the utility of Type Ia SNe as distance indicators and instead focus on the physics of an exploding white dwarf. As previously noted, PTF11kly was a very special supernova. In addition to exploding in a very nearby galaxy, PTF detected this SN just a few hours after it exploded, corresponding to the earliest detection of a Type Ia SN at the time. As Mansi highlighted in her talk, early detections of SNe can reveal a great deal about the progenitor systems. Here, we will look at how the PTF light curve of PTF11kly constrains the exploding white dwarf.
Theorists hate magnitudes, and prefer to work in the "natural" units of luminosity. Thus, to compare the observed PTF light curve to models we need to convert from magnitude to luminosity.
Problem A1 Convert the PTF11kly magnitude measurements to luminosity in units of $\mathrm{erg} \; \mathrm{s}^{-1}$, and store the results in an array called L
. Assume no bolometric correction from the $g$-band. Assume the PTF $g$-band has a central frequency of $6.284 \times 10^{14} \mathrm{Hz}$.
Hint - recall that AB magnitudes have a standard zeropoint:
$$m_{AB} = -2.5 \log_{10} \frac{f_\nu}{3631 \; \mathrm{Jy}}.$$
In [ ]:
# complete
# complete
# complete
L = # complete
Problem A2 Plot the luminosity, including upper limits, of PTF 11kly.
In [ ]:
plt.figure(figsize = (10,8))
plt.errorbar( # complete
plt.yscale( # complete
plt.xlim( # limit the plot to the few relevant upper limits prior to explosion
plt.xlabel( # complete
plt.ylabel( # complete
To compare models to the PTF 11kly light curve, we need to determine the exact time at which the SN exploded. Fortunately, the early luminosity evolution of Type Ia SNe has been shown to be parabolic:
$L(t) \propto (t - t_0)^2,$
where $L$ is the luminosity, $t$ is the time, and $t_0$ is the time of explosion.
Problem B1 Fit a parabolic function to the early (i.e., only fit the first few days) light curve to determine the time at which PTF11kly exploded.
In [ ]:
# complete
# complete
# complete
print('PTF11kly exploded at MJD = {:.3f}'.format( # complete
Now that we have determined the precise time of explosion, we can compare the luminosity of PTF11kly to theoretical models of shock breakout, which will constrain the radius of the progenitor. For example, Rabinak et al. (2011) found that the early luminosity of Type Ia SNe can be described as:
$$L(t) = 1.2 \times 10^{40} \frac{R_{10}E_{51}^{0.85}}{M_c^{0.69}\kappa_{0.2}^{0.85}f_p^{0.16}} t_d^{-0.31} \; \mathrm{erg} \; \mathrm{s}^{-1},$$where $R_{10}$ is the progenitor radius $R_p/10^{10} \; \mathrm{cm}$, $E_{51}$ is the explosion energy in units of $10^{51}\; \mathrm{erg}$, $M_c$ is the progenitor mass in units of $M_{\mathrm{ch}}$, and $\kappa_{0.2}$ is the opacity $\kappa/0.2 \; \mathrm{cm}^2 \; \mathrm{g}^{-1}$, and $f_p$ is the form factor.
Problem B2 Assuming $E_{51} = M_c = \kappa_{0.2} = 1$, and $f_p = 0.05$, plot the theoretical light curves for exploding white dwarfs with radii of $R_\mathrm{WD} = 0.01, 0.1, 1.0 \; R_\odot$.
In [ ]:
# complete
# complete
R_p = np.array( # complete
# complete
# complete
plt.figure(figsize = (10,8))
for r in R_p:
plt.plot( # complete
plt.legend( # complete
plt.yscale('log')
plt.xlim(0.0, 3.5)
plt.ylim( # complete
Rabinak et al. show that $L(t)$ is directly proportional to the progenitor radius, as should be confirmed by the above plot. From these curves, we can constrain the radius of PTF11kly.
Problem B3 Plot the early ($t \le 4$) light curve of PTF11kly on the above plot, along with the $t^2$-fireball model fit to show the comparison of the actual explosion to the models.
In [ ]:
# overplot these results on the theoretical models
plt.errorbar( # complete
plt.plot( # complete
Problem B4 Based on the above plot - what constraints can you place on the progenitor radius? Using the initial detection of PTF11kly, what is the maximum size of the progenitor?
In [ ]:
R = # complete
print('The radius of the progenitor is <= {:.3f} Rsun'.format( # complete
And thus, PTF11kly provides direct evidence that (at least one) Type Ia SN come from progenitors that are significantly smaller than the Sun! Something to chew on - is this proof that SNe Ia come from white dwarfs? Hint - can you think of any other astrophysical objects that are allowed within the constraints above?
If you finish early, work on the following problem, or continue working on this as homework for this evening.
The challenge problem is going to focus on improving the use of the PTF SNe for measuring cosmological distances. Typically, the goodness of an individual method is reported as the scatter (in mag) about the best fit Hubble line.
Challenge Problem 1 Assuming $H_0 = 73.24 \pm 1.74$, plot the Hubble expansion curve on a plot showing distance modulus, $\mu$, against redshift, $z$. Overplot the distance modulus and redshift of the SNe in the Firth et al. sample. Do you notice any trends?
In a separate plot, show the residuals relative to $H_0 = 73.24$, and calculate the scatter (rms in mag) of your method relative to this baseline. Modern SN light curve fitters produce a scatter of $\sim{0.14} \;\mathrm{mag}$. Given all that you now know - do you consider our method used to derive $H_0$ good?
In [ ]:
# complete
# complete
plt.figure()
plt.plot( # complete
# complete
# complete
plt.figure()
plt.plot( # complete
The most famous method for reducing the scatter is the Ia Hubble residuals involves correcting for the fact that more luminous SNe have slower evolving, or broader, light curves. This is known as the Phillips Relation, where the Type Ia absolute magnitude is correlated with the rate of decline (Phillips 1993). In addition to this, several other corrections have been searched for, with few proving to be as useful as the Phillips Relation.
Challenge Problem 2 Can you reduce the scatter in the Hubble residuals from Problem 1? There are several different corrections you could try in addition to something like $\Delta{m_{15}}$...
In [ ]: