Constraining Type Iax Supernova Progenitor Systems with Stellar Population Aging

TL;DR

This study uses stellar population analysis near Type Iax supernovae to constrain their delay times, favoring white dwarf progenitors over Wolf-Rayet stars based on local star ages.

Contribution

It provides the first empirical delay-time distribution for SNe Iax, constraining progenitor models through stellar population aging analysis.

Findings

Median delay time of approximately 63 million years.

Strong evidence favoring white dwarf progenitors over Wolf-Rayet stars.

Delay-time distribution constrains progenitor system characteristics.

Abstract

Type Iax supernovae (SNe~Iax) are the most common class of peculiar SNe. While they are thought to be thermonuclear white-dwarf (WD) SNe, SNe~Iax are observationally similar to, but distinct from SNe~Ia. Unlike SNe~Ia, where roughly 30\% occur in early-type galaxies, only one SN~Iax has been discovered in an early-type galaxy, suggesting a relatively short delay time and a distinct progenitor system. Furthermore, one SN~Iax progenitor system has been detected in pre-explosion images with its properties consistent with either of two models: a short-lived (<100 Myr) progenitor system consisting of a WD primary and a He-star companion, or a singular Wolf-Rayet progenitor star. Using deep \textit{Hubble Space Telescope} images of nine nearby SN~Iax host galaxies, we measure the properties of stars within 200 pc of the SN position. The ages of local stars, some of which formed with the SN progenitor system, can constrain the time between star formation and SN, known as the delay time. We compare the local stellar properties to synthetic photometry of single-stellar populations, fitting to a range of possible delay times for each SN. With this sample, we uniquely constrain the delay-time distribution for SNe~Iax, with a median and $1-\sigma$ confidence interval delay time of $63_{- 15}^{+ 58} \times 10^{6}$ years. The measured delay-time distribution provides an excellent constraint on the progenitor system for the class, indicating a preference for a WD progenitor system over a Wolf-Rayet progenitor star.

Figures

Tables

Table 1. Table 1: The list of observations used in this study. The data can be found here: http://dx.doi.org/10.17909/t9-qr61-xb59

SN	Instrument	Filters	Number of Exposures	Total Exposure	Date of Observation
Name			per Filter	Lengths (sec)	(UTC)
2008A	ACS/WFC	F555W, F625W, F775W	6, 6, 4	3750, 3530, 2484	08/18/2009
"	WFC3/IR	F110W	12	8336	08/18/2009
2008ge	ACS/WFC	F435W, F555W, F625W, F814W	2	1168, 768, 844, 1244	10/26/2012
2008ha	ACS/WFC	F435W, F555W, F625W, F814W	2	9068, 764, 840, 12058	01/02/2013, 12/30/2015
2010ae	ACS/WFC	F435W, F555W, F625W, F814W	2	1402, 1002, 1078, 1478	05/23/2014
2010el	ACS/WFC	F435W, F555W, F625W, F814W	2	1168, 768, 844, 1244	05/23/2013
"	WFC3/UVIS	F275W, F336W, F438W, F555W, F814W	3, 3, 3, 5, 5	2382, 1119, 965, 1853, 1769	09/2013, 07/2015, 08/2017
2011ay	WFC3/UVIS	F555W, F814W	2	780, 780	01/20/2018
2012Z	ACS/WFC	F435W, F555W, F814W	8, 10, 10	9624, 12642, 12868	09/2006
2014ck	WFC3/UVIS	F625W, F814W	2	510, 680	02/22/2013
2014dt	WFC3/UVIS	F275W, F438W	2, 20	858, 400	02/28/2017, 11/18/2014

Equations

$$\mathrm{IGoF}=\frac{\sqrt{{\sum }_i{\mathrm{PhysDist}}_i\cdot min({\sum }_k(m_{ik}-m_{Iso.,k}{)}^2)}}{{\sum }_i{\mathrm{PhysDist}}_i},$$

$$\mathrm{PhysDist}$$

$$\mathrm{Normal}(\mu =0,\sigma =.65\frac{\mathrm{km}}{\mathrm{s}}\cdot \mathrm{Age}),$$

$$\displaystyle p(t)=\left\{\begin{array}[]{cc}0&\hskip 14.22636ptt<a\ \text{Myrs},\\ t^{-b}&\hskip 14.22636ptt>=a\ \text{Myrs}.\\ \end{array}\right.$$

$$\mathrm{Like}(\mathrm{Data},D,R|Iso.)$$

$$=\prod _i{L}_i,$$

$$L_i$$

$$L_{ij}$$

$$p(Iso.|Data)=\iint$$

$$p(Iso.|Data,D,R)$$

$$p(Iso.)$$

$$p(Age|Data)=\int p(Iso.|Data)p(M)dM,$$

Full text

††pubyear: 2018

1 Introduction

Type Iax supernovae (SNe Iax) are a class of peculiar SNe which share some characteristics with Type Ia SNe, but appear to be physically distinct (Foley et al., 2013). SNe Iax are also known as SN 2002cx-like SNe, after the prototypical object of the class (Li et al., 2003; Jha et al., 2006; Jha, 2017). These SNe are principally characterized by their low luminosities (as compared to SNe Ia), low photospheric velocities, a lack of a secondary maximum in near-IR bands (Li et al., 2003), and their unique late-time spectra (e.g., Jha et al., 2006; Foley et al., 2016). In contrast to SNe Ia, SNe Iax reach their peak brightness in optical bands before they do in the near-IR.

SNe Iax are the most common type of peculiar SN, occurring at roughly 30% the rate of SNe Ia (Foley et al., 2013; Miller et al., 2017). This, along with their observational similarities with SNe Ia, means that investigating the physical properties of these supernovae is extremely helpful to understanding SNe Ia. The comparably low photospheric velocities of SNe Iax makes line identification easier, allowing for precise measurements of the explosion in these events, which can help us to understand the physics involved in both SNe Iax and Ia. This will be crucial as we look to improve our knowledge of SNe Ia in order to continue to perform precision cosmology. Despite decades of study of SNe Ia, the identity of the binary companion to the SN Ia progenitor is still unknown. By first learning about the progenitor system for SNe Iax, we may be able to shed some light on this long-standing problem.

As compared to typical SNe Ia, SNe Iax are also inferred to have much lower explosion energies (Branch et al., 2004), and lower ejecta masses (Foley et al., 2013) (though this has been debated for some members of the class, e.g., Sahu et al. 2008, Stritzinger et al. 2015, Yamanaka et al. 2015). Additionally, SNe Iax tend to have a layered structure in their ejecta, similar to SNe Ia. In contrast to SNe Ia, however, there is evidence of Ni mixing in at least two SNe Iax, further from the center than in standard SNe Ia (Jha et al., 2006; Phillips et al., 2007; Stritzinger et al., 2015). Finally, at least six SNe Iax show cobalt in their near-IR spectra (Kromer et al., 2013; Stritzinger et al., 2014, 2015; Tomasella et al., 2016). Taken together, these properties suggest that SNe Iax may be the result of the partial deflagration of a white dwarf (WD) (e.g. Foley et al., 2010b; Kromer et al., 2013; Fink et al., 2014; Magee et al., 2016).

SNe Iax host galaxies are also markedly different from SNe Ia host galaxies. While the SNe are nearly always found in late-type galaxies, there is one example (SN 2008ge) of a SN Iax in an S0 galaxy (Foley et al., 2010a). Lyman et al. (2013) and Lyman et al. (2018) find evidence of star forming regions near the sites of most SNe Iax, suggesting a young progenitor age at the time of explosion for these events. Additionally, SNe Iax seem to preferentially occur in metal-poor regions of their host galaxies (Lyman et al., 2018). This is quite different from the metallicity of the explosion sites for SNe Ia, instead matching quite closely with the sites of SNe Ib, Ic, II, and IIb. Though the metallicities of the explosion sites for SNe Iax and core-collapse SNe are similar, they appear quite different from those of low-redshift long gamma-ray bursts (LGRBs; Levesque et al. 2010). This is particularly interesting, as the most popular core collapse model for SNe Iax is the “fallback SN” model (Valenti et al., 2009; Moriya et al., 2010), which are theorized to occur alongside LGRBs (Heger et al., 2003; Della Valle et al., 2006). Thus the divergent metallicity distributions of SNe Iax and LGRBs suggest that SNe Iax may not be “fallback” core-collapse SNe.

The physical origin of SNe Iax is still uncertain, but the leading models suggest that they are the result of binary interactions between a carbon-oxygen WD and a He-star companion (Foley et al., 2013; Jha, 2017). These interactions seem to result in an incomplete deflagration (sub-sonic nuclear burning, Phillips et al. 2007) of the WD which may not completely unbind the star (e.g., Jordan et al., 2012; Kromer et al., 2013; Fink et al., 2014; Long et al., 2014; Magee et al., 2016). The low luminosities, velocities, and ejecta masses measured from these events together give strong indirect evidence for a WD progenitor system (Foley et al., 2009, 2010b). Additionally, the Ni mixing observed in some SNe Iax ejecta can be most easily explained by a turbulent deflagration in a WD (Röpke, 2005). Despite the mounting evidence for this model, there remain several unresolved issues. Specifically, the incomplete deflagration model struggles to reproduce the observed lower luminosity SNe Iax, such as SNe 2008ha and 2010ae (Kromer et al., 2013; Fink et al., 2014) It has also been argued (Kromer et al., 2015) that SN 2008ha spectra are indicative of poorly mixed ejecta, in contrast to brighter SNe Iax. Well-mixed ejecta are a hallmark of partial deflagration models, making SN 2008ha particularly confusing. Additionally, at least two SNe Iax (SNe 2004cs and 2007J) have He lines in their spectra, which is difficult to reproduce with the WD and He-star binary model (Fink et al., 2014).

For the WD He-star binary progenitor model for SNe Iax, it is expected that the circumstellar environment will be enriched by pre-SN mass loss, either from winds from the He donor star, or from non-conservative mass transfer. Once the WD explodes, the interaction between the blast wave and the circumstellar material should power X-ray and radio emission (Chevalier & Fransson, 2006; Immler et al., 2006; Russell & Immler, 2012; Margutti et al., 2012, 2014). Though X-ray emission has not yet been detected from any SNe Iax, Liu et al. (2015) use X-ray upper limits from seven SNe Iax (Russell & Immler, 2012) to constrain their pre-SN mass loss rates. In comparing theoretical pre-SN mass loss rates to the observed X-ray upper limits, they find broad agreement across a variety of models and observations.

The WD He-star binary model for SNe Iax has also been tested with pre-explosion Hubble Space Telescope (HST) imaging of the locations of SN 2008ge (Foley et al., 2010a), SN 2008ha (Foley et al., 2014), SN 2012Z (McCully et al., 2014a), and SN 2014dt (Foley et al., 2015). These works use precision astrometry to align pre-explosion images with ground-based images of the SN. The authors then perform photometry on the pre-explosion images at the site of the explosion to probe the progenitor system just before explosion. In doing so, they find an upper bound on the progenitor mass for SN 2008ge (assuming a massive star progenitor), find an upper bound on the age of the SN 2008ha progenitor system, and resolve the likely progenitor system of SN 2012Z. Foley et al. (2015) find no progenitor system for SN 2014dt down to quite deep limits ($3-\sigma$ limits of $M_{F438W}>-5.0$ mag and $M_{F814W}>-5.9$ mag), making the interpretation of SN 2014dt as a core-collapse event less likely. While these limits rule out most Wolf-Rayet star models, there are some Wolf-Rayet star models which remain below the detection limits. In all four cases, the results show consistency with the WD and He-star model, though the large photometric errors involved prevent them from ruling out other models with high certainty.

The likelihood of observing a WD progenitor in pre-explosion imaging of SN Iax explosion sites is quite low, due to the inherent low luminosity of WDs as compared to their mass donor companion. However, analyzing the host environments of SNe (in pre or post-explosion imaging) has been quite successful in extracting information on the progenitors of Type Ib, Type Ic, and Type II SNe (e.g., Maíz-Apellániz et al., 2004; Badenes et al., 2009; Kuncarayakti et al., 2013a, b). One such method uses resolved stellar populations near the SNe, fit to theoretical isochrones (lines of constant age on a color-magnitude diagram) to measure the age of the SN progenitor (e.g., Gogarten et al., 2009; Murphy et al., 2011; Jennings et al., 2012; Williams et al., 2014; Jennings et al., 2014). In this method, it is assumed that the resolved stellar population and the progenitor itself formed nearly simultaneously, in a single burst of star formation. This method of aging stellar populations to measure time between star formation and explosion – known as the delay-time – for SNe has only been effective for SNe with relatively short delay times, as typical open cluster velocity dispersions are large enough that stellar positions are no longer correlated after a few hundred million years (Bastian & Goodwin, 2006; Lada, 2010). Every model for SNe Iax predicts a delay-time distribution (DTD) with non-zero probability down to below 100 Myrs, so it is reasonable to apply this method to measure SNe Iax delay times.

In the WD progenitor model for SNe Iax, explosions occur at the Chandrasekhar mass, and the quickest binary channel for a carbon/oxygen (C/O) WD is to accrete helium from a He star companion (Hachisu et al., 1999; Postnov & Yungelson, 2014). The stable mass transfer rate can be high for helium accretion, and Claeys et al. (2014) show that this channel dominates the thermonuclear SN rate between 40 Myr and 200 Myr (above which traditional SNe Ia dominate). This has been observed in several binary population synthesis studies (Ruiter et al., 2009; Ruiter et al., 2011; Wang et al., 2009b, a; Meng & Yang, 2010; Piersanti et al., 2014; Liu et al., 2015). Liu et al. (2010) present a model (originally intended to explain a different kind of system) that begins with a 7 M⊙ + 4 M⊙ binary system that undergoes two phases of mass transfer and common envelope evolution, before resulting in a 1 M⊙ C/O WD + 2 M⊙ He star. As the He star evolves, it again fills its Roche lobe and begins stable mass transfer onto the white dwarf that could lead to the SN Iax. The high accretion rate involved in this process mean that total delay time is dominated by the stellar evolution timescale of the secondary, giving short expected delay times for these models.

In this work, we employ Bayesian Monte Carlo methods to fit isochrones to the stellar populations around SNe Iax, in order to accurately measure the DTD of this class of SNe. With a sample of nine SNe Iax, we constrain the ages of these systems at the time of explosion.

Li et al. (2018) examine the location of SNe Iax in their host galaxies, and find that they tend to occur at much larger projected radii as compared to SNe Ia. The authors perform a Kolmogorov-Smirnov test to compare the fractional host galaxy fluxes at the explosion sites for SNe Iax, SNe Ia, and SNe Ib/Ic. They find strong evidence that SNe Ia and SNe Iax are drawn from different populations, while also finding that they cannot reject the hypothesis that SNe Iax and SNe Ib/Ic are drawn from the same fractional host galaxy flux distribution. Taken together with the metallicity information, this suggests that the delay-time for SNe Iax might be much more similar to those of core-collapse SNe than the delay-times measured for SNe Ia.

This paper is structured in the following way. In Section 2, we describe the observations used in this study, and the methods to extract photometry of the stellar populations near each SN. In Section 3, we detail the Bayesian Monte Carlo methods used to probabilistically determine the ages of these stellar populations. We discuss the particulars of each SN individually in Section 4, describing the priors used in the fit and the resulting posterior distributions. In Section 5, we discuss the overall properties of the measured DTDs for this sample of SNe Iax, and describe how our results compare to theoretical models for this class of SN.

2 Observations and Data Reduction

2.1 Sample Selection

The distance to a host galaxy is the key parameter which restricts our sample. Because the methods we will employ depend on analyzing a population of well resolved stars in the direct vicinity of each SN under consideration, only relatively nearby targets are suitable. Because of this limitation, our 2013 HST program (GO-12999) observed SNe 2008ge (Pignata et al., 2008), 2008ha (Foley et al., 2009), 2010ae (Pignata et al., 2010), and 2010el (Monard, 2010). At the time our HST program was executed, these were the four closest known SNe Iax. They were selected as a sample to demonstrate that stellar population fitting could be a reliable method to probe the delay-time of SNe Iax. These four objects were observed with the Advanced Camera for Surveys Wide Field Channel (ACS/WFC), with the F435W, F555W, F625W, and F814W filters. These bands were chosen because they provide good colors for observing broad spectral types, and their spectral range enables one to correct for extinction due to dust in the local environment.

In the time since 2013, we have searched the HST archive for other SNe Iax with high quality data (wide spectral range and long exposures) of their stellar neighborhoods. This search added SNe 2008A (GO-11590: Jha, S.), 2011ay (GO-15166: Filippenko, A.), 2012Z (GO-10497: Riess, A.; GO-10711: Noll, K.; GO-10802: Riess, A.; and GO-13757: Jha, S.), 2014ck (GO-13029: Filippenko, A.), and 2014dt (GO-13683: Van Dyk, S.; GO-14779: Graham, M.) to our sample. Two of these objects have ACS/WFC imaging, with SN 2012Z imaged in the F435W, F555W, and F814W bands, and SN 2008A imaged in the F555W, F625W, and F775W bands. SN 2011ay and SN 2014ck each have Wide Field Camera 3 (WFC3) imaging, in the F555W and F814W bands, and the F625W and F814W bands, respectively. Additional imaging of SNe 2010el was added in this search (GO-13364: Calzetti, D.; GO-13816: Bentz, M.; GO-14668: Filippenko, A.; and GO-15133: Erwin, P.). Finally, our 2015 follow-up HST program for SN 2008ha (GO-14244: Foley, R.) provided additional imaging in the F435W and F814W filters.

Our complete sample of SNe Iax contains all observed SNe Iax within 35 Mpc, for which HST host galaxy imaging is deep enough to find strong constraints on progenitor age. In addition, two more distant SNe (2008A and 2011ay) are included in the sample, thanks to the abundance of HST imaging of their host galaxies. SNe 2008A, 2008ge, 2008ha, 2010ae, 2010el, 2011ay, 2012Z, 2014ck, 2014dt compose our sample. The complete list of data can be found in Table 1. Of these, previous works have found constraints on the progenitors of SNe 2008ge, 2008ha, and 2012Z. These constraints were achieved through pre-explosion HST imaging of the sites of SNe 2008ge (Foley et al., 2010a) and 2012Z (McCully et al., 2014a), and post-explosion imaging of the site of SN 2008ha (Foley et al., 2014), with claims of detections of either the remnant, or the probable donor star for both SN 2008ha and SN 2012Z.

Note that for host galaxy distances which are measured using redshift, we assume $H_{0}=73.24\pm 1.74$ km Mpc*-1* s*-1* (Riess et al., 2016). The distances to each object, as well as the methods with which these distance were measured, can be found in Section 4. The Additionally, Milky Way reddening is assumed to follow the extinction maps of Schlafly & Finkbeiner (2011), using $R_{v}=3.1$ (Cardelli et al., 1989). This reddening is applied to all simulated data to match the observations.

2.2 Data Reduction

To build catalogs from the HST observations, we use a custom pipeline written primarily in Python111https://github.com/cmccully/snhst. The pipeline initially registers the astrometry to a ground-based image which has a wider field of view than the HST images. The World Coordinate System (WCS) from the ground-based image is considered to be the global astrometric solution. This stage typically produces a precision of $\sim 0.05\arcsec$ which corresponds to roughly one HST pixel for ACS.

After solving for the global WCS, we refine the registration between individual exposures, starting from the flat-fielded frames from the MAST archive222https://archive.stsci.edu/. For ACS and WFC3, we use the FLC frames that have been corrected for charge-transfer inefficiency (CTI) using the pixel-based method (Anderson & Bedin, 2010). Cosmic rays are rejected from the individual frames using Astro-SCRAPPY333https://github.com/astropy/astroscrappy before the registration process to alleviate false positives in the catalog matches. After doing a coarse registration between HST visits by hand, the pipeline uses TweakShifts from Drizzlepac444https://drizzlepac.stsci.edu (STSCI Development Team, 2012) to refine the offsets between frames.

Once the registration is completed, we combine the exposures using Astrodrizzle adopting standard values for the parameters${}^{\ref{foot1}}$. We include the cosmic-ray step in Astrodrizzle to do final cosmic-ray rejection. For cases that we have more than 4 individual exposures in the same filter for a target, we subsample the pixel grid and decrease the drizzle pixel coverage fraction to 0.8.

To build the final catalogs, we run DOLPHOT, a modified version of HSTPhot (Dolphin, 2000), using the drizzled image as a coordinate reference. DOLPHOT runs on the individual flat fielded frames and stacks the photometry to produce the final catalog, including all point sources across the entirety of each image. We again use the CTI-corrected FLC frames for ACS and WFC3/UVIS. For WFPC2, we use the CTI correction built into DOLPHOT. DOLPHOT includes its own image registration stage which generally produces a scatter of $\sim 0.01\aas@@fstack{\prime\prime}$ Finally, we inject artificial stars using DOLPHOT to estimate the brightness limit of our images.

2.3 Star Selection

Using the DOLPHOT output, we apply the recommended cuts in sharpness (<0.3), roundness (<1), and crowding (<0.1), along with a minimum signal-to-noise ratio of 3.5. Though this removes most non-stellar detections, there remain a number of clusters and bright clumps of gas which are selected as stars by DOLPHOT. As such, for each region surrounding a SN, every detected source within the region is checked by eye to ensure that it has a point spread function indicative of a single stellar source. This allows us to remove gas clumps from our catalogues.

As the human eye struggles to distinguish distant clusters from stars in the galaxy under examination, another method is used to make this distinction. Aperture photometry is performed using IRAF in order to determine a concentration parameter for each source (Chandar et al., 2010). The concentration parameter gives the difference in measured F555W (or a similar band) magnitude between when a 3 pixel aperture is used, and when a .5 pixel aperture is used. By separating on this concentration parameter, as detailed in Chandar et al. (2010), we remove the extended sources from our sample, leaving us with only sources that have a high likelihood of being stars.

3 Isochrone Fitting Method

In order to determine the age of a supernova progenitor in in post-explosion host galaxy images, we assume a star formation history (SFH) characterized by a dominant star formation event for the region in which the progenitor system formed. We also assume a typical scale within which we can reasonably expect nearly every star examined to have a shared SFH. Past work has indicated that this radius is anywhere from 50 pc to 200 pc (Bastian & Goodwin, 2006; Eldridge et al., 2011; Williams et al., 2014; Foley et al., 2014; Maund & Ramirez-Ruiz, 2016; Maund, 2017). To account for this uncertainty, we consider all stars within a 200 pc projected on-sky radius of the SN position, weighting the stars according to their probability of being associated with the SN.

To estimate the probability of association as a function of on-sky distance, we build a probability distribution as follows. We first assume a flat initial distribution of stars within 100 pc of the SN, to represent the cluster at the time of formation. We then assume a velocity dispersion of $0.65$ km/s (Geller et al., 2009), multiplied by the age of the cluster (in this case the delay-time being tested), to model the cluster spreading out over time. To “apply” this spreading effect to the initial flat distribution, we convolve the two distributions. This convolution of the two distributions gives us a rough probability of association with the SN as a function of projected physical distance. This probability function is used to weight the stars within 200 pc, with the higher probability stars receiving a higher weight in the fitting scheme. We fit a single stellar population model to each of these populations.

Geller et al. (2009) quotes a velocity dispersion for clusters of $0.65\pm 0.10$ km/s. To account for this uncertainty, we run our complete analysis using dispersions of $0.55$ km/s, $0.65$ km/s, and $0.75$ km/s. Each of these choices ultimately leaves the final PDF in delay time largely unchanged. For this reason, we use the median value of $0.65$ km/s in all the analyses that follow.

Maund (2017) uses a similar technique to age the stellar population around the positions of 12 Type IIP SNe. In that paper, however, they fit to multiple stellar populations for each SN, allowing for multiple bursts of star formation in the SFH. In this study, we choose to only fit a single stellar population around each SN, as the errors in magnitude and color for our detected stars are large enough that fitting to multiple populations would result in overfitting.

3.1 Color Magnitude Diagrams

To date the regions around the SNe, we use the magnitudes and colors derived using DOLPHOT to place each star on a color magnitude diagram (CMD). We then overplot the MIST synthetic photometry isochrones (Paxton et al., 2011; Paxton et al., 2013, 2015; Dotter, 2016; Choi et al., 2016), corrected for distance, metallicity, and extinction. We test isochrones in the age range $10^{6.5}$ to $10^{8.5}$ years. For metallicity and extinction, we test the $2$-$\sigma$ range from our priors. One sigma errors in distance and internal (host galaxy) reddening are propagated through to color and magnitude, added in the appropriate way and displayed in CMDs as shaded regions.

3.2 Isochrones

In order to fit the synthetic photometry isochrones to the detected stars accurately and with well understood errors, a hierarchical Bayesian framework is used to create a new statistic to measure goodness of fit555https://github.com/TTakaro/Type-Iax-HST. Monte Carlo methods are then used to translate this statistic into probability distributions. In analogy to a Chi-squared fit, this statistic (which we will hereafter call IGoF for “Isochrone Goodness of Fit") measures the minimized difference in magnitude in each filter between the stars detected near the SN and the isochrone being tested, with the magnitude difference in each filter then summed in quadrature. Summing over each detected star, IGoF gives a total distance from the stars to the isochrone in n-dimensional magnitude space, where n is the number of filters. As described in the previous section, each star is weighted in IGoF according to its projected physical distance from the SN, and thus its likelihood of forming coincident with the SN progenitor system. This statistic was selected in part because of its qualitative similarity to a Chi-squared fit, allowing for relatively simple analysis.

The equation for IGoF is shown below.

[TABLE]

[TABLE]

where $i$ specifies a given star, $k$ specifies a given filter, $m_{ik}$ is the magnitude of a given star in a filter, $m_{\mathrm{Iso},jk}$ is the magnitude drawn from the isochrone, adjusted for distance, milky way reddening, and host galaxy reddening. $\mathrm{PhysDist}_{i}$ is the weight given to the star $i$ according to its physical distance from the SN. As detailed in section 3, this weight is the convolution of a flat distribution, and a normal distribution, evaluated at the on-sky distance measured between star $i$ and the SN. Because a cluster of stars disperse as the cluster ages, this weighting function will have a different shape at different ages.

To translate the measured value of IGoF for an isochrone into a relative probability that the stars in the data match that isochrone, careful forward modeling is required. A set of artificial stars equal in number to the detected stars are generated from each isochrone, drawn from normal distributions in flux space using the characteristic flux error from the detected stars in the host galaxy image. These stars are given radial distances on the sky generated from the radial distribution used to weight stars in the IGoF. An IGoF value is then measured for the set of artificial stars. This process is repeated 5000 times for each isochrone, until the histogram of IGoF values converges to a probability distribution. This probability density function (PDF) gives a probability of measuring a value of IGoF, given the chosen age, metallicity, and host galaxy reddening associated with the isochrone. Using this fit PDF, the values of IGoF measured from the data are translated into values of relative probability for each isochrone. These relative probabilities are then normalized to determine a relative probability distribution in age, metallicity and host galaxy extinction for each SN. We then marginalize over metallicity and host galaxy extinction in order to extract the one dimensional age distributions for each object.

Further details of our method, including the likelihood function used in the isochrone fitting scheme are shown in Appendix A.

3.3 False Star Tests and Photometric Completeness

To determine the completeness of our photometry, we perform false star recovery tests using the false star tool in DOLPHOT. For host galaxy image, we insert 50,000 stars with magnitudes covering the whole range of measured magnitudes, and with x-y positions covering the entirety of the region around the SN. Using the results of the false star runs in DOLPHOT and binning with 0.1 mag resolution, we calculate a recovery fraction as a function of magnitude for each host galaxy image. We then apply this recovery fraction to the artificial stars generated while calculating IGoF in the process mentioned above.

In several of the imaged SN host galaxies, no stars are detected within 200 pc of the SN position. In order to get an upper bound on the age of the progenitor star for these SNe, we use the recovery fraction calculated using false star testing. From this recovery fraction function, we determine a limiting magnitude of detection in the image, by requiring that 90% of inserted stars be recovered for each magnitude bin above our limiting magnitude. We then simulate the effects of a $50M_{\odot}$ open cluster (Lada 2010 find $50M_{\odot}$ to be the most common open cluster mass) at the position of the SN under consideration, drawing stars from a Kroupa Initial Mass Function (IMF) (Kroupa, 2001) until we reach the total cluster mass. If any of these drawn stars have brightnesses above our limiting magnitude given the isochrone we are considering, we rule out the associated age, metallicity, and reddening combination for this SN. By performing this analysis for each isochrone, we establish a lower bound on the age of each SN progenitor system for which we detect no nearby stars (SNe 2008A, 2008ge, and 2011ay).

4 Analysis of Individual Objects

The objects in our sample and the associated data used for this study are listed in Table 1.

4.1 SN 2008A

SN 2008A was discovered in NGC 624 (Nakano et al., 2008), at an estimated Tully-Fisher distance of $51.5\pm 11.0$ Mpc (Theureau et al., 2007). For this object, we assume no host galaxy extinction, as the SN is located on the fringes of it’s host galaxy (Lyman et al., 2018). The metallicity for this object is difficult to measure (Lyman et al., 2018), so we use solar metallicity as our prior (Asplund et al., 2009). For additional analysis of SN 2008A, see e.g., Milne et al. (2010), Hicken et al. (2012), McCully et al. (2014b), and Foley et al. (2016).

In our photometric analysis, we find no stars within 200 parsecs of the SN position. As such, we use the false stars method of DOLPHOT to set a lower limit on the delay time for this object. This analysis measures a minimum delay time for SN 2008A of $\rm Age_{08A}\geq 6$ Myr.

4.2 SN 2008ge

SN 2008ge was discovered in NGC 1527 (Pignata et al., 2008), at an estimated Tully-Fisher distance of $17.37\pm 0.96$ Mpc (Tully et al., 2013). We assume no host galaxy extinction for this object (Foley et al., 2010a). For this object, we also use a solar metallicity prior (Jorgensen, 1997).

In our photometric analysis, we find no stars within 200 parsecs of the SN position. As with SN 2008A, we use the false stars method of DOLPHOT to set a lower limit on the delay time for this object. From this analysis, we find $\rm Age_{08ge}\geq 18$ Myr.

4.3 SN 2008ha

SN 2008ha was discovered in UGC 12682 (Puckett et al., 2008), which has an redshift-derived distance of $21.3\pm 1.5$ Mpc (Foley et al., 2014). Lyman et al. (2018) report a metallicity of [Fe/H] $=-0.78\pm 0.09$ in the region directly around the SN. Additionally, Foley et al. (2014) report no indication of host extinction in the region around the supernova. For more information on SN 2008ha, see e.g., Valenti et al. (2009), Foley et al. (2009), Pumo et al. (2009), Fryer et al. (2009), Foley et al. (2010b), Stritzinger et al. (2014), and Foley et al. (2016). Our analysis uses 16 stars in the vicinity of SN 2008ha to perform the isochrone fitting.

Using the measurement from Lyman et al. (2018) as a prior, metallicity is allowed to float between [Fe/H] = $-0.50$ and $-1.00$. Our fitting algorithm returns a metallicity posterior which is slightly skewed towards solar metallicity, as compared to the prior. Our fits show no preference for any host galaxy reddening for the stars around the SN, consistent with the literature (Foley et al., 2014). Marginalizing over metallicity, we measure the delay time for SN 2008ha of $52_{-9}^{+13}$ Myr. The full probability distribution in delay time is shown in Figure 1.

4.4 SN 2010ae

SN 2010ae was discovered in ESO 162-17 (Pignata et al., 2010), which is determined using the Tully-Fisher relation to be $11.09\pm 1.02$ Mpc away (Tully et al., 2013). Lyman et al. (2018) report a metallicity of [Fe/H] $=-0.43\pm 0.06$ in the region directly around the SN. A high (though uncertain) upper limit on host extinction of $E(B-V)_{\textrm{host}}=0.50\pm 0.42$ mag is reported in the region around the supernova (Stritzinger et al., 2014). Additional information on SN 2010ae can be found in e.g., Foley (2013) and Foley et al. (2016). We find 11 stars in the vicinity of SN 2010ae, which we use for our isochrone fitting.

We use Gaussian priors for both metallicity and host-galaxy extinction. The posterior distribution in extinction is skewed towards low extinction values, with a strong preference for $E(B-V)_{\textrm{host}}\leq 0.36$ mag. The posterior distribution in metallicity on the other hand is completely consistent with the prior distribution. Marginalizing over both metallicity and host-galaxy extinction, we find a delay time of $117_{-29}^{+19}$ Myr for SN 2010ae. The full probability distribution for the age of the progenitor is shown in Figure 2. The measured probability distribution in age has two peaks, one at 9 Myrs, the other at 110 Myrs. This is likely indicative of two separate stellar populations at the explosion site of SN 2010ae, either of which may have been associated with the SN. Providing further evidence for two separate stellar population in the area is the preferred host galaxy extinction values for the isochrone fit. The younger peak favors high values for extinction, while the older peak favors lower values, indicating that the younger population is in the background as compared to the older, as its light passed through more dust.

The probability distribution in delay time that we measure for this object is particularly illuminating, as it demonstrates that our method is able to recover the short delay times that we would expect from a WR progenitor scenario. The young stars around SN 2010ae are likely a stellar population that is unrelated to the SN progenitor, but our detection of this population indicates that our analysis pipeline is working correctly.

4.5 SN 2010el

SN 2010el was discovered in NGC 1566 (Monard, 2010), at a Tully-Fisher distance of $5.63\pm 1.12$ Mpc (Tully et al., 2013). Lyman et al. (2018) report a metallicity of [Fe/H] $=0.16\pm 0.07$ at the site of the SN. No precise measurement is available of the host galaxy extinction in the region around the SN, though a measurement of the extinction of the galaxy at large has been measured to be $E(B-V)_{\textrm{host}}=0.205$ mag (Gouliermis et al., 2017). We use 46 stars in the vicinity of SN 2010el to perform our isochrone fitting.

Though we analyse both ACS/WFC data and WFC3/UVIS data for this object, DOLPHOT is better able to identify stars in the ACS/WFC3 data. As such, our ACS/WFC data is far more discernining than the WFC3/UVIS data which we analyze. For this reason, the PDF in age for this object is determined from the ACS/WFC data. All of the analysis for SN 2010el that follows is based on this ACS/WFC data.

We use the metallicity estimate for the SN as a prior, while assuming a flat prior of $0\leq E(B-V)_{\textrm{host}}\leq 0.205$ mag for the extinction. Our posterior in metallicity is consistent with the prior, showing a slight preference towards higher metallicity. Our posterior in host galaxy extinction shows a strong preference for non-zero extinction, with the preferred value of $E(B-V)_{\textrm{host}}=0.205$ mag. When we marginalize over metallicity and host-galaxy extinction, we find a delay time of $53_{-6}^{+5}$ Myr for SN 2010el. The age probability distribution is shown in Figure 3.

4.6 SN 2011ay

SN 2011ay was discovered in NGC 2315 (Blanchard et al., 2011), to which we measure a distance of $87.4\pm 6.4$ Mpc, derived from the galaxy’s redshift (Miller & Owen, 2001). Using the measured metallicity gradient for NGC 2315, we infer a metallicity of [O/H] = $-0.19\pm 0.20$ (Lyman et al., 2018). The SN has a relatively large offset from its host galaxy, with the SN occurring in a rather empty environment. As such, we assume no host extinction as our prior. For more detailed analysis of SN 2011ay, see White et al. (2015), Szalai et al. (2015), Foley et al. (2016), and Barna et al. (2017).

As with SN 2008A and SN 2008ge, in our photometry of NFC 2315, we find no stars within 200 parsecs of the SN location. We perform the same false star tests as with these two objects, in the hope of getting a lower limit on the delay time. However, because of the great distance to NGC 2315, our method cannot rule out any of the isochrones considered in this study.

4.7 SN 2012Z

SN 2012Z was discovered in NGC 1309 (Cenko et al., 2012), which is determined to be $31.92\pm 0.88$ Mpc away (Riess et al., 2016), using the Cepheid distance method. Lyman et al. (2018) report a metallicity of [Fe/H] $=-0.43\pm 0.08$ at the site of the SN, in agreement with metallicity gradient that they measure. Host galaxy reddening is estimated to be $E(B-V)_{\textrm{host}}=0.07+-0.03$ mag using Na I and K I doublets (Stritzinger et al., 2015). More information on SN 2012Z can be found in e.g., McCully et al. (2014a), Yamanaka et al. (2015), and Foley et al. (2016). We use 10 stars in the vicinity of SN 2012Z to perform our isochrone fitting.

Using these measurements as Gaussian priors, we fit for both metallicity and host galaxy extinction. In each case, the calculated posterior is consistent with the prior, showing a slight preference for metallicity and host reddening on the upper end of the ranges given in the literature. Marginalizing over metallicity and extinction, we measure a delay time of $61_{-11}^{+13}$ Myr for SN 2012Z. The full probability distribution in age is shown in Figure 4.

4.8 SN 2014ck

SN 2014ck was discovered in UGC 12182 (Masi et al., 2014), which is determined to be $24.32\pm 1.69$ Mpc away (Tomasella et al., 2016), based on its redshift. Taddia et al. (2015) report a metallicity of [Fe/H] $=-0.51\pm 0.26$ at the site of SN 2006fp, which is also located in UGC 12182, at roughly the same projected distance from the galaxy nucleus. The host galaxy extinction is estimated to be $E(B-V)_{\textrm{host}}\leq 0.05$ mag, determined using the Na I absorption line (Tomasella et al., 2016). We use 20 stars in the vicinity of SN 2014ck to perform our isochrone fitting.

These measurements are used as priors, with the metallicity and host galaxy extinction then fit and the posterior in metallicity and extinction measured. Our posterior in metallicity is peaked around the same value as the prior, though with noticeably thinner tails, indicating that our data show a strong preference for [Fe/H] $=-.51$. Our data also indicate no preference for higher extinction, showing consistency with $E(B-V)_{\textrm{host}}=0$ mag. Marginalizing over both metallicity and extinction, we measure a delay time of $113_{-25}^{+21}$ Myr for SN 2014ck. The full age probability distribution is shown in Figure 5.

4.9 SN 2014dt

SN 2014dt was discovered in Messier 61 (Ochner et al., 2014) at a distance of $19.3\pm 0.6$ Mpc (Rodríguez et al., 2014) from our Milky Way galaxy. Lyman et al. (2018) report a metallicity of [Fe/H] $=0.16\pm 0.07$ at the explosion site of SN 2014dt, in rough agreement with the metallicity of [Fe/H] $=-0.01$ reported in Foley et al. (2015). Foley et al. (2015) also report no indications of host-galaxy reddening at the site of the SN. We adopt these as priors for our fit.

Our method finds no likely stars within 200 pc of the SN. Performing the false star analysis described in Section 3.3, we are unable to confidently rule out any ages that we test.

5 Discussion

5.1 Verification of Fitting Method

As mentioned earlier, the probable progenitor star for SN 2012Z has previously been detected (McCully et al., 2014a), while a progenitor search for SN 2008ha has yielded strong limits on the delay time for this object (Foley et al., 2014). In Foley et al. (2014), the authors also fit Padova isochrones to the stars directly surrounding the probable remnant of SN 2008ha, to measure an approximate delay-time for the SN. They find a 1-$\sigma$ confidence interval on the age of the system of $55_{-10}^{+13}$ Myr. Using our isochrone fitting method, we find a 1-$\sigma$ confidence interval of $52_{-9}^{+13}$ Myr, in excellent agreement with the literature. In McCully et al. (2014a), the authors similarly fit isochrones by eye to the stars surrounding the possible progenitor of SN 2012Z, deriving a best fitting age range of 10 - 42 Myrs. With our analysis, we find a 1-$\sigma$ confidence internal of $61.0_{-11}^{+13}$ Myr. This disagreement is only at the $\sim 1\sigma$ level and likely comes from the fact that McCully et al. (2014a) used several sources in close proximity to the SN in their fit, which our analysis rejects as being too extended to confidently classify as stars. Our fit instead relies on sources which are further from the SN, but which are more likely to truly be stars. This new preference for higher ages now strongly suggests that SN 2012Z is not consistent with a WR progenitor system.

5.2 Comparisons Between Objects

We check for correlations between delay time and peak luminosity, and between delay time and metallicity (as (Liu et al., 2015) suggest). Our peak luminosities are taken from the Open Supernova Catalog (Guillochon et al., 2017), while our metallicities are the priors used for our fitting (citations listed above). For this check, we perform Pearson r-tests to measure the strength of a correlation (or anti-correlation) between our measured delay time and either metallicity or peak luminosity. We find r values very close to zero in each case, indicating no statistically significant correlations.

Though SN 2008ge is unique in its host galaxy characteristics, it does not appear to stand out in this study. Instead, we find that the lower limit which we measure from SN 2008ge is completely consistent with the measured DTDs of all but one of the SN in our study. The DTD for SN 2010ae is double peaked, showing two stellar populations. The older of these two populations is consistent with the lower bound from SN 2008ge, though the younger population is not. Because of this broad agreement between SN 2008ge and the sample at large, for the purposes of this study SN 2008ge appears to belong to the class of “standard” SNe Iax. SN 2014ck has also been called an outlier among SNe Iax due to its observed properties (Tomasella et al., 2016). However our analysis does not find that its delay time in significantly different from the delay times of the rest of our sample.

5.3 SN 2014dt Progenitor Mass Limit

The false star tests that we perform for SN 2014dt give us a $3-\sigma$ magnitude limit of 24 mag in the F438W filter. Though this is not sufficient to rule out any isochrones that we test, it is sufficient to rule out single stars more massive than $19.5M_{\odot}$ at the distance of the host galaxy of SN 2014dt. Using MESA single stellar evolution models of the appropriate metallicity, this roughly corresponds to a lower limit on the delay time for this particular object of $10$ Myrs. Because this limit is not derived using the same methods as the rest of the study, we do not include the limit in the rest of the analysis.

5.4 Empirical Measurements

If we assume that all SNe Iax that we consider in this study share a progenitor channel, we can make strong statistical statements about the delay-time distribution of SNe Iax at large. The probability density function (PDF) for each SN, and the summed PDF (unnormalized for clarity in the figure) are shown in Figure 6. Taking only the SNe for which a full DTD can be measured (ignoring lower limits), we find a median delay time and $1-\sigma$ confidence interval of $63_{-15}^{+58}$ Myrs. The data show a strong preference for a delay time peaked near 55 Myrs, with the PDFs for three SNe peaking within 6 Myrs of 55 Myrs. SN 2014ck and SN 2010ae on the other hand, represent the long tail of the distribution to higher ages, with the medians of their DTDs at 113 Myr and 117 Myr respectively.

If we then draw groups of 5 SNe (the number of SNe Iax in the summed PDF) from our summed PDF and fit power-law decay models to the result in an MCMC fashion, we find the posterior distribution in Figure 7. This is done using emcee (Foreman-Mackey et al., 2013) with a flat prior on the cutoff age (0 to 100 Myr) and an inverse-gamma prior with a peak at -1 for the decay exponent. The resulting form of this power-law decay is then described by Equation 2.

[TABLE]

The distributions of the cutoff age (labeled $a$) and the decay power (labeled $b$) are shown in Figure 7. The 1-sigma confidence intervals for each parameter are $40.83_{-13.38}^{+7.85}$ Myrs and $2.52_{-0.81}^{+1.18}$ respectively.

As seen in Figure 7, the distribution in cut-off age ($a$) is double peaked. This double peaked distribution is due to the double peaked measured DTD for SN 2010ae. The less probable, younger peak in $a$ is due to the fit occasionally drawing an age from the younger peak of the SN 2010ae DTD, and trying to fit this alongside the older ages from the other SNe. This could just be the result of a 2nd stellar population forming near the explosion site of SN 2010ae by chance, or it could be indicative of a heterogeneity in SNe formation channels. Repeating this study with a larger sample of SNe Iax – for instance, once many more are discovered by the Large Synoptic Survey Telescope (LSST) – would allow for one to answer this question one way or another.

5.5 Model Comparisons

In Valenti et al. (2009), Moriya et al. (2010), and Pumo et al. (2009), the authors suggest a core-collapse origin for at least some SNe Iax. These models predict a maximum delay time for SNe Iax of 6 to 10 Myrs (Fryer, 1999; Heger et al., 2003). If we ignore the events with no nearby stars detected, we find just a $2\%$ probability of at least one SN in our sample having a delay time of 10 Myr or less, strongly suggesting a long-lived progenitor channel for SNe Iax. Using only the analysis of SN 2008ge, we find a $3-\sigma$ minimum age of 18 Myrs, giving us strong evidence that SNe Iax are not produced by short-lived progenitors.

The WD He-star binary model for the progenitor systems of SNe Iax instead predicts a delay-time of roughly 30 Myrs and up, with an expected peak in probability at the lower bound delay-time, and a long tail up to much higher ages (Wang et al., 2014; Liu et al., 2015). The power law decay which we fit to the summed PDF above is a good parameterization of this expected DTD. The best-fit lower limit on delay-time for our data is $40.83$ Myrs, in excellent agreement with the model predictions.

6 Conclusions

We have developed a new method to fit simulated photometry of stellar populations generated from MIST isochrones to broadband photometric measurements, weighted by distance from a cluster center. We then use this method to fit the stellar neighborhoods of nine SNe Iax. We employ Bayesian methods to generate probability distributions for the delay times for these nine objects. The probability distributions are shown in Figure 6. We find a 68% confidence interval for the delay time of our sample of $63_{-15}^{+58}$ Myrs. When fitting a power law to the overall probability distribution, we find a 68% confidence interval on the lower bound on the delay time of $40.83_{-13.38}^{+7.85}$ Myrs.

Taken together, this sample of SNe Iax provides the most precise measurement of the DTD for SNe Iax to date. The data strongly disfavor a fallback core-collapse origin for SNe Iax, instead showing consistency with a WD-He star binary progenitor model.

Acknowledgements

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with programs GO-12999 and GO-14244.

The authors would like to thank A. Skemer for his helpful advice in the design of the Isochrone Goodness of Fit statistic used in this paper. We also thank J. Maund and C. Kilpatrick for their helpful advice on using DOLPHOT. We thank R. Hounsell, and D. Kasen and L. Bildsten for their contributions to one of the two HST proposals that led to this paper. We additionally thank A. Wasserman for his help with the final statistical analyses of the DTD, and J. Schwab for a variety of helpful discussions. Finally, we thank R. Murray-Clay, without whose advice and patience this work would not have been possible.

The UCSC team is supported in part by HST programs GO-12999 and GO-14244, NSF grant AST-1518052, the Gordon & Betty Moore Foundation, the Heising-Simons Foundation, and by a fellowships from the David and Lucile Packard Foundation to R.J.F. This SN Iax research at Rutgers University is supported by NSF award 161545 and HST programs GO-11133, GO-11590, GO-12913, and GO-13757. M. Stritzinger is supported by a research grant (13261) from the VILLUM FONDEN.

Appendix A Additional Information on Isochrone Fitting

The likelihood function which we use in our isochrone fitting scheme is shown below.

[TABLE]

where D is distance, R is host galaxy reddening, $i$ specifies a given star, $j$ specifies a given stellar mass in the isochrone, $k$ specifies a given filter, $m_{ik}$ is the magnitude of a given star in a filter, $\sigma_{ik}$ is the corresponding uncertainty, and $m_{\mathrm{Iso},jk}$ is the magnitude from an isochrone in a given filter for a given mass. $m_{\mathrm{Iso},jk}$ is the magnitude drawn from the isochrone, adjusted for distance, Milky Way reddening, and host galaxy reddening.

The posterior distribution for an isochrone given a stellar population is then:

[TABLE]

Here $p(\mathrm{Age})$ and $p(\mathrm{Metallicity})$ are the prior distributions in age and metallicity. Our prior in age is a flat distribution between $10^{6.5}$ and $10^{8.5}$ years. Our metallicity priors are described for each SN in section 4. This posterior distribution for an isochrone is the relative probability that the stellar population being considered has a specific age and metallicity, given the stellar data.

Our measured DTD for each SN is this posterior distribution, marginalized over metallicity.

[TABLE]

where M is metallicity.

Appendix B Example of Isochrone Fitting

We will now walk through the process of fitting a set of isochrones to a data set of stars around a supernova. For this example, we will use SN 2008ha. As mentioned in Section 2, we run DOLPHOT on a set of HST images of the host-galaxy of SN 2008ha. We apply the recommended list of cuts on the data output in order to get a list of point sources in the image. This list is restricted to include only those stars which are within 200 pc projected on-sky distance from the SN, assuming a distance to the host-galaxy of $21.3\pm 1.5$ Mpc. Clusters are removed using the method detailed in Section 2.3 in order to get our final list of stars which we will use for isochrone fitting.

The library of isochrones used in this study are logarithmically spaced in age, using 0.05 dex increments from $10^{6.5}$ to $10^{8.5}$ years. These isochrones are drawn using the median, $1-\sigma$, $2-\sigma$ metallicities of our prior, and adjusted for the internal reddening of the host-galaxy. We then construct a distribution to measure probability of a star forming with the SN progenitor, as explained at the beginning of Section 3. Now for each star in our data, we apply the weighting distribution evaluated at the distance of the star. We then find the minimum of the square root of the distance from the star to each isochrone, summed across all of the filters in the data. This process is completed for each star in the data, summing the values measured for each star in quadrature. These values (different for each age, metallicity, and extinction) are our measured Isochrone Goodness of Fit (IGoF) values.

To convert the measured IGoF values to probability distributions in age, metallicity and host-galaxy reddening, we perform careful forward modelling. Stars are drawn from a Kroupa IMF, using the isochrone under consideration to convert the mass to a magnitude in each filter in our data, keeping each star if its associated magnitudes are above our observational limit. Once there are the same number of synthetic stars as real stars in the data, we apply a distance error to all of the synthetic stars together by drawing from a normal distribution with parameters determined by our distance prior. These stars are shown in the second panel of Figure 8. Magnitude errors are applied by drawing from a normal distribution for each star and filter, using the measured magnitude uncertainties as a function of magnitude which are derived from the data. These stars appear in the third panel of Figure 8. We now measure a value of IGoF, and repeat this process 5000 times for each isochrone, in order to get a probability distribution in IGoF for each isochrone. This PDF is used to relate our measured value of IGoF for the isochrone to a relative probability associated with the age, metallicity, and reddening of the isochrone.

The fourth panel of Figure 8 shows the real stars detected in the vicinity of SN 2008ha, to compare to the artificial stars which we draw from the isochrone. Note that the real data has higher dispersion than the artificial data because the real data likely features stars which were not born in the same star forming event as the SN progenitor. These stars would have been drawn from isochrones of a different age, and hence would look like a bad fit for the isochrone we plot, whereas every artificial star is drawn from the same isochrone that is plotted.

In order to further verify our method, we generate false stars from our isochrone, applying characteristic errors to the stars and our observational sensitivities, and employ our method to recover the age and metallicity of the isochrone. The results of this effort are shown in Figure 9, Figure 10, and Figure 11. These figures have injected metallicities that correspond to the median from our prior ([Fe/H] $=-0.78$), while the injected ages are $10^{7.45}$, $10^{7}$, and $10^{7.8}$ years. These tests are successful, as we recover the correct age as the approximate peak in each distribution. Note that each of these tests is performed with 15 stars, in order to match the data for SN 2008ha. For SNe around which we detect more stars (as in SN 2010el, for example), our method recovers the correct age with even higher precision.