Think Global, Act Local: The Influence of Environment Age and Host Mass on Type Ia Supernova Light Curves

TL;DR

This study investigates how the age of the host galaxy environment and local factors influence Type Ia supernova light curves, revealing significant correlations that impact cosmological measurements like the Hubble constant.

Contribution

The paper introduces an improved Bayesian method to estimate host environment ages and demonstrates their significant effect on supernova luminosity residuals, refining understanding of supernova standardization.

Findings

A 0.114 mag step in Hubble residuals at ~8 Gyr host age.

A new principal component, PC₁, correlates strongly with residuals (4.7σ).

At most a 1.3% impact on H₀ estimation.

Abstract

The reliability of Type Ia supernovae (SNIa) may be limited by the imprint of their galactic origins. To investigate the connection between supernovae and their host characteristics, we developed an improved method to estimate the stellar population age of the host as well as the local environment around the site of the supernova. We use a Bayesian method to estimate the star formation history and mass weighted age of a supernova's environment by matching observed spectral energy distributions to a synthesized stellar population. Applying this age estimator to both the photometrically and spectroscopically classified Sloan Digital Sky Survey II supernovae (N=103) we find a $0.114 \pm 0.039~{\rm mag}$ `step' in the average Hubble residual at a stellar age of $\sim 8~\text{Gyr}$; it is nearly twice the size of the currently popular mass step. We then apply a principal component analysis on the SALT2 parameters, host stellar mass, and local environment age. We find that a new parameter, PC$_1$, consisting of a linear combination of stretch, host stellar mass, and local age, shows a very significant ($4.7\sigma$) correlation with Hubble residuals. There is a much broader range of PC$_1$ values found in the Hubble flow sample when compared with the Cepheid calibration galaxies. These samples have mildly statistically different average PC$_1$ values, at $\sim 2.5\sigma$, resulting in at most a 1.3% reduction in the evaluation of H$_0$. Despite accounting for the highly significant trend in SNIa Hubble residuals, there remains a 9% discrepancy between the most recent precision estimates of H$_0$ using SNIa and the CMB.

Figures

Tables

Table 1. Table 1: Local environment data for Campbell et al. ( 2013 ) SNIa

SDSS ID	redshift	uncert.	u	g	r	i	z	${\sigma }_u$	${\sigma }_g$	${\sigma }_r$	${\sigma }_i$	${\sigma }_z$	HR	uncert.
		$\times {10}^{-5}$	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]
762	0.19138	2.4	24.65	23.82	22.95	22.61	22.18	0.79	0.09	0.04	0.04	0.13	0.15	0.08
1032	0.12975	3.3	24.92	24.74	23.73	23.32	22.87	0.61	0.16	0.09	0.11	0.26	-0.15	0.12
1371	0.11934	2.6	23.22	21.52	20.42	20.00	19.62	0.16	0.01	0.00	0.006	0.01	-0.13	0.06
1794	0.14191	6.3	23.76	23.77	23.09	22.89	22.81	0.43	0.11	0.08	0.06	0.24	0.27	0.08
2372	0.18046	2.2	24.81	22.85	21.84	21.40	21.02	0.87	0.03	0.01	0.01	0.03	-0.12	0.06

Table 2. Table 2: Data used for validation with Gupta et al. ( 2011 ) ages

SDSS ID	redshift	u	g	r	i	z	${\sigma }_u$	${\sigma }_g$	${\sigma }_r$	${\sigma }_i$	${\sigma }_z$
		[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]
1166	0.3820	22.54	21.83	20.04	19.38	19.03	0.48	0.11	0.04	0.03	0.08
1580	0.1830	24.99	20.41	19.21	18.72	18.28	1.35	0.03	0.02	0.02	0.04
2165	0.2880	22.82	23.35	22.04	22.22	21.28	0.38	0.24	0.12	0.20	0.36
2422	0.2650	23.64	22.86	22.00	21.95	22.01	0.90	0.23	0.17	0.24	0.76
2789	0.2905	22.01	20.92	19.42	18.84	18.39	0.35	0.05	0.02	0.02	0.06

Table 3. Table 3: Global host data for Campbell et al. ( 2013 ) SNIa

SDSS ID	u	g	r	i	z	${\sigma }_u$	${\sigma }_g$	${\sigma }_r$	${\sigma }_i$	${\sigma }_z$
	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]	[mag]
762	20.34	18.50	17.46	17.01	16.70	0.13	0.01	0.01	0.01	0.02
1032	21.49	19.40	18.30	17.83	17.47	0.19	0.01	0.01	0.01	0.02
1371	20.60	18.62	17.55	17.10	16.68	0.08	0.01	0.01	0.01	0.01
1794	22.37	20.76	20.37	20.10	20.12	0.45	0.05	0.06	0.08	0.29
2372	21.79	20.50	19.53	19.02	18.59	0.24	0.03	0.02	0.02	0.05

Table 4. Table 4: SEDs and redshifts for nearby galaxies

			Local SED					Global SED
Host Galaxy	SNIa	redshift	u	g	r	i	z	u	g	r	i	z
M101	2011fe	0.000804	14.73	13.28	12.87	12.65	12.71	11.63	10.14	9.56	9.24	9.02
NGC 1015	2009ig	0.008797	14.54	12.93	12.19	11.78	11.53	18.49	16.91	16.19	15.74	15.54
NGC 1309	2002fk	0.007138	15.10	13.96	13.46	13.18	13.01	12.93	11.97	11.54	11.33	11.19
NGC 3021	1995al	0.00535	15.24	14.05	13.44	13.12	12.86	13.76	12.61	12.02	11.71	11.44
NGC 3370	1994ae	0.004276	17.08	15.97	15.50	15.23	15.06	13.63	12.58	12.06	11.79	11.56
NGC 3447	2012ht	0.003556	17.53	16.33	15.88	15.65	15.55	14.78	13.82	13.41	13.23	13.07
NGC 3972	2011by	0.002799	16.85	15.69	15.09	14.71	14.60	14.07	12.86	12.24	11.88	11.66
NGC 3982	1998aq	0.00371	15.28	14.16	13.72	13.43	13.27	12.95	11.92	11.45	11.17	10.98
NGC 4424	2012cg	0.00162	13.93	12.90	12.35	11.88	11.85	12.97	11.87	11.28	10.78	10.77
NGC 4536	1981B	0.00603	16.51	15.30	14.77	14.50	14.43	12.60	11.32	10.64	10.28	10.03
NGC 4639	1990N	0.00364	18.04	16.91	16.46	16.14	16.12	13.40	12.11	11.46	11.12	10.89
NGC 5584	2007af	0.005525	17.01	15.94	15.47	15.16	15.26	13.39	12.41	11.92	11.59	11.64
NGC 7250	2013dy	0.0039	15.35	14.63	14.17	14.08	13.88	14.27	13.48	13.04	12.93	12.72
UGC 9391	2003du	0.00649	17.80	16.79	16.39	16.21	16.06	16.05	15.15	14.78	14.62	14.42

Table 5. Table 5: SFH parameters used in the self consistency test

ID	$\log \left(z/z_{☉}\right)$	${\tau }_2$	$\tau$	$t_0$	$t_i$	$\varphi$	age
			[Gyr-1]	[Gyr]	[Gyr]	[rad]	[Gyr]
1	-0.5	0.1	0.5	1.5	9.0	-0.785	10.68
2	-0.5	0.1	0.5	1.5	9.0	1.504	1.41
3	-0.5	0.1	7.0	3.0	10.0	1.504	1.75
4	-0.5	0.1	7.0	3.0	13.0	0.0	4.28
5	-1.5	0.1	0.5	1.5	9.0	-0.785	10.68
6	-0.5	0.8	7.0	3.0	10.0	1.504	1.75
7	-0.5	0.1	0.5	1.5	6.0	1.504	2.40
8	-0.5	0.1	0.1	8.0	12.0	1.52	0.437

Table 6. Table 6: SEDs for the self consistency test

ID	u	g	r	i	z	$⟨A⟩$
	[mag]	[mag]	[mag]	[mag]	[mag]	[Gyr]
1	20.36	18.76	17.99	17.67	17.39	$8.5\pm 1.5$
2	20.31	18.74	17.98	17.66	17.39	$7.7\pm 1.5$
3	16.15	15.43	15.40	15.19	15.21	$1.4\pm 0.5$
4	17.65	16.74	16.49	16.26	16.16	$4.2\pm 1.0$
5	19.69	18.29	17.70	17.45	17.29	$7.2\pm 1.6$
6	17.66	16.58	16.25	16.01	15.86	$2.7\pm 0.8$
7	17.62	16.80	16.57	16.34	16.26	$4.5\pm 1.4$
8	19.72	18.37	17.88	17.68	17.56	$4.4\pm 1.1$

Table 7. Table 7: Results of median age, stellar mass, and PC 1

SN	local age	${\sigma }_a$	global age	${\sigma }_a$	$\log \left(\text{M}/{\text{M}}_{\odot }\right)$	PC1
762	5.1	2.7	8.0	3.9	11.1	0.07
1032	5.8	2.6	8.6	4.3	10.5	-1.78
1371	8.9	4.3	8.1	1.3	10.7	-0.82
1794	4.0	1.4	2.0	0.9	8.8	2.20
2372	5.9	1.5	4.8	2.9	10.2	-0.07

Table 8. Table 8: Normalization parameters applied before PCA

	$\mu$	$\sigma$
$x_1$	-0.177	1.015
$c$	0.0100	0.0829
$\log \left(M/M_{\odot }\right)$	10.15 dex	0.69 dex
age	5.22 Gyr	2.11 Gyr

Table 9. Table 9: PCA coefficients using local and global ages

	$x_1$	$c$	$\log \left(\text{M}/{\text{M}}_{\odot }\right)$	ageaathe median age of the local environment posterior	% variance	$x_1$	$c$	$\log \left(\text{M}/{\text{M}}_{\odot }\right)$	agebbthe median global age posterior	% variance
$P{C}_1$	0.56	-0.10	-0.54	-0.63	44	0.45	-0.13	-0.60	-0.64	47
$P{C}_2$	-0.16	0.96	-0.21	-0.12	25	-0.21	0.94	-0.26	-0.10	25
$P{C}_3$	-0.65	-0.26	-0.71	0.07	18	-0.85	-0.32	-0.39	-0.19	19
$P{C}_4$	0.49	0.09	-0.40	0.77	11	-0.16	0.07	0.65	-0.74	8

Table 10. Table 10:

	$x_1$	$c$	mass	age
$\mu$	0.54	-0.09	-0.52	-0.61
$\sigma$	0.08	0.17	0.11	0.06
min	-0.29	-0.95	-0.90	-0.89
max	0.85	0.93	0.31	0.06

Table 11. Table 11: Ages, stellar mass, and SALT2 values from the H 0 subscript H 0 \text{H}_{0} SNIa calibration sample

SNIa	host	$x_1$	${\sigma }_{x_1}$	$c$	${\sigma }_c$	citation	local age	${\sigma }_a$	global age	${\sigma }_a$	$\log \left(\text{M}/{\text{M}}_{\odot }\right)$	PC1,local	PC1,global
1981B	NGC 4536	-0.32	0.14	0.030	0.010	J07	6.5	2.4	3.0	2.4	10.2	-0.49	0.61
1990N	NGC 4639	0.63	0.04	0.014	0.004	J07	6.5	1.4	4.8	3.0	10.2	0.06	0.55
1994ae	NGC 3370	0.32	0.10	-0.065	0.033	J07	5.4	1.3	5.0	1.9	10.0	0.46	0.62
1995al	NGC 3021	0.71	0.08	0.051	0.006	J07	6.4	1.7	6.0	1.5	10.2	0.06	0.16
1998aq	NGC 3982	-0.40	0.07	-0.086	0.007	J07	5.2	1.4	6.1	1.4	10.0	0.10	0.01
2002fk	NGC 1309	0.22	0.04	-0.101	0.003	S12aaLight curve data supplied by the Open Supernova Catalog (Guillochon et al., 2016).	6.6	2.6	5.2	1.4	10.1	-0.06	0.44
2003du	UGC 9391	0.30	0.04	-0.100	0.004	J07	4.6	1.3	4.4	1.4	9.0	1.50	1.73
2007af	NGC 5584	-0.45	0.02	0.053	0.004	H12	6.2	1.6	6.2	2.0	9.8	-0.26	-0.13
2009ig	NGC 1015	1.76	0.15	-0.058	0.013	H12	8.8	4.4	7.9	2.1	10.3	0.01	0.20
2011by	NGC 3972	0.02	0.13	0.012	0.014	B14aaLight curve data supplied by the Open Supernova Catalog (Guillochon et al., 2016).	4.5	3.0	4.1	2.6	9.8	0.60	0.78
2011fe	M 101	-0.21	0.07	-0.066	0.021	P13bbCitation is for the SALT2 parameters.	3.1	1.0	4.9	1.2	9.9	0.88	0.46
2012cg	NGC 4424	0.45	0.04	0.080	0.020	V18bbCitation is for the SALT2 parameters.	3.5	2.3	2.8	2.5	10.0	0.88	1.07
2012ht	NGC 3447	-1.25	0.05	-0.080	0.030	V18bbCitation is for the SALT2 parameters.	4.0	1.3	4.4	1.4	9.2	0.63	0.79
2013dy	NGC 7250	0.70	0.04	0.089	0.025	V18bbCitation is for the SALT2 parameters.	4.6	1.2	4.2	1.4	9.2	1.33	1.51

Equations

$$(Z{e}^{-Z}{)}^{\text{pmetals}}$$

$$Z\equiv \frac{z}{z_{\text{sun}}10^{\log (z/z_{\text{sun}})}}$$

$${\tau }_{\lambda }(t)=\{\begin{array}{ll}{\tau }_1(\lambda /5500\text{A˚}{)}^{-0.7}& t\le 10^7\text{yr}\\ {\tau }_2(\lambda /5500\text{A˚}{)}^{-0.7}& t>10^7\text{yr}\end{array}$$

$${\Psi }_{\star }(t)\propto \{\begin{array}{ll}(t-t_0)e^{-(t-t_0)/\tau }& t_0\le t\le t_i\\ {\Psi }_{\star }(t_i)+m_{\text{sf}}\mathcal{R}(t-t_i)& t_{\text{i}}<t\le \mathcal{A}(z)\\ 0& \text{else}\end{array}$$

$$⟨A{⟩}_{\text{mass}}=\mathcal{A}(z)-t_0-\frac{{\int }_{t_0}^{\mathcal{A}(z)}(t-t_0){\Psi }_{\star }(t)dt}{{\int }_{t_0}^{\mathcal{A}(z)}{\Psi }_{\star }(t)dt}$$

$$⟨A{⟩}_{\text{mass}}=\mathcal{A}(z)-t_0-\frac{{\int }_0^{\mathcal{A}(z)-t_0}t{\Psi }_{\star }^{\prime }(t)dt}{{\int }_0^{\mathcal{A}(z)-t_0}{\Psi }_{\star }^{\prime }(t)dt}$$

$$P(\theta |D)=\frac{P(D|\theta )P(\theta )}{P(D)}$$

$$P(\theta |D)\propto P(\theta )\times \mathcal{L}(D|\theta ).$$

$$\ln (\mathcal{L})\propto \sum _i\frac{(m_i-(m_{\text{FSPS},i}+\delta ){)}^2}{{\sigma }_i^2}+\ln \left(2\pi {\sigma }_i^2\right)$$

$$2.5$$

$$0.1$$

$$-1.520838$$

$$0.5$$

$$-45.0$$

$$x_1^{\prime }=(x_1-{\mu }_{x_1})/{\sigma }_{x_1}.$$

$${\mathrm{PC}}_1=0.56x_1^{\prime }-0.10c^{\prime }-0.54m^{\prime }-0.63a^{\prime },$$

$$\mathrm{HR}=0.051\text{mag}\times {\mathrm{PC}}_1-0.012\text{mag},$$

$$\mu =m_B-M_B+(\alpha -{\alpha }^{\prime })x_1-\beta c+\gamma m^{\prime }+\gamma a^{\prime }$$

$$10^{\Delta M/5}-1.$$

Full text

distance scale, galaxies: distances and redshifts, galaxies: general, galaxies: photometry, galaxies: stellar content, supernovae: general

††journal: ApJ††software: astropy (Astropy Collaboration, 2013), astroquery (10.6084/m9.figshare.805208), corner.py (Foreman-Mackey, 2016), emcee (Foreman-Mackey et al., 2012), FSPS (Conroy et al., 2009; Conroy & Gunn, 2010), kcorrect (Blanton & Roweis, 2007), Matplotlib (Hunter, 2007), Numpy (van der Walt et al., 2011), Pandas (McKinney, 2010), Python-FSPS (Foreman-Mackey et al., 2014), scikit-learn (Pedregosa et al., 2011), SciPy (Jones et al., 2001), Seaborn (10.5281/zenodo.883859), sfdmap (), SIMBAD (Wenger et al., 2000), sncosmo (10.5281/zenodo.592747)

1 Introduction

For decades, astronomers have been developing methods to better understand the variation in peak luminosities of Type Ia supernovae (SNIa) and improve their use as precision distance indicators. In 1993, Phillips identified a relationship between peak magnitude and the rate of fading in the light curves of SNIa. A connection between supernova color and peak luminosity was also shown to improve distance estimates of SNIa (Riess et al., 1996; Tripp & Branch, 1999; Phillips et al., 1999). SNIa quickly became useful cosmological probes, measuring the expansion rate of the universe (Hamuy et al., 1995; Riess et al., 1995), the density of matter (Garnavich et al., 1998; Perlmutter et al., 1998), and the accelerated expansion of the universe (Riess et al., 1998; Perlmutter et al., 1999).

The principle energy source that powers a SNIa light curve derives from the decay chain of ${}^{56}\text{Ni}$ that is synthesized in the runaway nuclear fusion at the start of the explosion (Arnett, 1982). The radioactive nickel yield appears to vary by a factor of eight over the extreme range of typical SNIa luminosities. The mass of radioactive elements helps regulate the rate of recombination in iron group elements, and this results in the famous “Phillips relation” between the light curve decline rate and luminosity (Kasen & Woosley, 2007).

The reason some SNIa synthesize nearly a solar mass of radioactive nickel while others are powered by a tenth of a solar mass remains uncertain. Models suggest that the speed of the fusion front moving through the white dwarf has a major influence on the radioactive yield. The transition between deflagration (subsonic fusion) and detonation (supersonic fusion) may vary from supernova to supernova and this could explain the diversity in their luminosities. The variation in nickel yield appears to be influenced by host properties as Hamuy et al. (1996) and Hamuy et al. (2000) noted that low-luminosity SNIa tend to occur in passive galaxies like large ellipticals. This observation has been confirmed and expanded in several subsequent studies (Gallagher et al., 2005; Sullivan et al., 2006) which clearly demonstrate that host galaxies provide an important clue to the progenitors and explosion mechanisms of SNIa.

Metal abundance and population age are global properties of galaxies that could conceivably have an impact on a supernova’s ${}^{56}\text{Ni}$ yield. Host properties that correlate with age or metallicity, such as galaxy mass, could also influence the character of a supernova explosion. For example, Timmes et al. (2003) proposed that the metal abundance during the main sequence stage could affect the neutron fraction in the resulting white dwarf stars. They predicted that high metal abundance populations will generate low radioactive yields and, thus, faint supernovae. Applying the galactic mass-metallicity relationship to this finding, SNIa from large galaxies would then be systematically fainter just as was seen in Hamuy et al. (2000). Bravo et al. (2010) derived a similar luminosity-metallicity relationship, but 3D calculations by Röpke & Hillebrandt (2004) suggest the effect is much smaller than originally envisioned. An observational test of this hypothesis by Gallagher et al. (2008) looked at SNIa in elliptical galaxies and found trends with both age and metallicity, although this was disputed by Howell et al. (2009).

A sensitive test of the environmental effects on SNIa is to study the scatter in a SNIa Hubble diagram after light curve shape and color corrections have been applied. This type of analysis is also important for constraining systematic errors in cosmological measurements. Hubble residuals are the difference between the luminosity distance and the expected distance from the best fit cosmology and are most powerful in the “Hubble flow” where peculiar velocities are small compared with the expansion velocity. Woosley et al. (2007) and Kasen et al. (2009) showed that explosion asymmetries, metallicity variations, kinetic energy variations, and other explosion parameters can produce a dispersion in the Phillips relation for a fixed radioactive yield. These relationships present the possibility of using Hubble residuals to probe supernova physics. Research over the past several years indicates that some host galaxy properties correlate with Hubble residuals (e.g. Sullivan et al., 2010; Lampeitl et al., 2010; Childress et al., 2013).

Surprisingly, the strongest correlation between Hubble residuals and galaxy properties appears to be with the host stellar mass. The effect is called the “mass step” because at $\sim 10^{10}~{}\text{M}_{\sun}$ there appears to be a jump in the average Hubble residual. This led to extensive work on understanding the physical properties such as population age and metallicity that could underlie the mass correlation (e.g. Gupta et al., 2011; Hayden et al., 2013; Moreno-Raya et al., 2016a, b). Childress et al. (2014) has proposed that the mass correlation is really an age effect amplified by galaxy “downsizing.” This research has been fruitful, but not definitive. Ongoing analyses of SNIa data sets continue to debate the significance of these effects (e.g. Graur et al., 2016; Jones et al., 2018, for LOSS and Foundation respectively).

Rigault et al. (2013) and Rigault et al. (2018) looked at star formation and specific star formation rates respectively at the sites of SNIa explosions by measuring $\text{H}_{\alpha}$ emission strength using spatially resolved spectroscopy. They found a very significant step in corrected SNIa peak luminosity as a function of the specific star formation rate at the location of the explosion. This research identified a set of SNIa with a 0.2 mag offset in Hubble residual that are generally found in regions of lower star formation rate. There is still not a consensus on the impact of local star formation rates on Hubble residual especially when these trends are measured using UV observations (Rigault et al., 2015; Jones et al., 2015).

A very recent study (Jones et al., 2018) on a large low-redshift data set (the Foundation sample, Foley et al., 2017) compared Hubble residual with host galaxy stellar mass, local environment mass (the region within a radius of $1.5~{}\text{kpc}$ of the SNIa), host galaxy rest frame $u-g$ color, and local environment rest frame $u-g$ color. The rest frame $u-g$ color is a simple way to estimate recent star formation, and therefore, a crude age estimator. They found that the local environment contained no new information compared to the global parameters.

In addition, there appears to be a tension between the $\text{H}_{0}$ estimated from the cosmic microwave background observations (Planck Collaboration, 2016; Planck Collaboration et al., 2018) and estimates based on SNIa. These precision measurements now disagree by $3.8\sigma$ (Riess et al., 2016, 2018) due to either new physics or a systematic bias in one of the measurements. The SNIa host environment is of particular importance to the precision measurements of the Hubble constant ($\text{H}_{0}$). The Cepheid variables used to calibrate the SNIa peak absolute magnitude are massive stars, so are observed only in star-forming galaxies. In contrast, supernovae discovered in unbiased searches are found in all types of galaxies. A mass step correction between the Cepheid calibrated hosts and the Hubble flow galaxies is currently applied to the data (at $\pm 0.03~{}\text{mag}$) and provides a relatively minor tweak to the value of $\text{H}_{0}$ derived from supernovae. So, at present, host environment is not a major contributor to the uncertainty budgets of cosmological measurements.

Here, we scrutinize the relation between Hubble residual and the age of the stellar population derived from host galaxy colors. We develop a technique to translate multi-band galaxy photometry into an estimate of the star formation history and finally an average age for the stellar population. In principle, colors should be a better indicator of SNIa progenitor age than $\text{H}_{\alpha}$. This is because after a single burst of star formation, $\text{H}_{\alpha}$ emitting HII regions are dissipating just as SNIa are beginning to explode. We apply our technique to both the global photometry of SNIa host galaxies and to the specific populations at the sites of the explosions. We then compare the local and global ages with SNIa characteristics, and other host properties to determine the parameters that have the largest impact on the measured Hubble residual. Finally, we investigate how our estimated ages may influence the current SNIa measurements of $\text{H}_{0}$ and its error budget.

2 Data

For our primary analysis, we use SNIa that were discovered with the SDSS-II Supernova survey (Sako et al., 2008).We selected both spectroscopic and photometric classified SNIa that passed cosmology cuts as described and released by Campbell et al. (2013). This cosmological data set, including their SALT2 (Guy et al., 2007, 2010) calibration parameters, arehttp://www.icg.port.ac.uk/stable/campbelh/SDSS_Photometric_SNe_Ia.fits (catalog available online). Hereafter we will refer to Campbell et al. (2013) as C13. For our analysis, a few additional cuts were applied. For quality, we limited our analysis to objects whose photometric uncertainty is less than $1.5~{}\text{mag}$. These cuts are dominated by low quality $u$-band magnitudes. In the end, the resulting $g$- and $i$-band maximum uncertainties are less than $0.3~{}\text{mag}$, and less than $0.16~{}\text{mag}$ for $r$-band. In addition we removed objects that had a Hubble residual greater than $0.7~{}\text{mag}$. Looking at Figure 16 of C13, these objects are likely core collapse supernovae that passed the color-magnitude cut.

The software developed for the data aggregation and analysis described in this paper is available at https://github.com/benjaminrose/mc-age.

2.1 Local Environment Photometry

The photometry of the environment around the site of the supernovae are taken from the “Scene Modeling” method described in Holtzman et al. (2008). The method incorporates all the images taken during the survey to build a photometric model at the location of the transient. The resulting photometry at the site is a convolution of the galaxy on the scale of the typical seeing. By applying this fixed angular aperture we always get the most compact local environment possible. Finally, in order to keep the environment truly “local,” a redshift cutoff was imposed, $z<0.2$. With the average seeing for SDSS being $1.4^{\prime\prime}$, the maximum extent of a galaxy’s local environment was $3~{}\text{kpc}$ in radius. The typical size of the projected aperture defining the environment at the supernova location was $1.5~{}\text{kpc}$.

This results in a final data set of 103 SNIa. A partial list of the final data set used in the local environment analysis is visible in Table 1.

2.2 Global Photometry

In addition to the photometry of the local environment, we also study the correlation between the supernova characteristics and the host properties as a whole (hereafter the global properties). From the global host properties we can compare our results directly to the work presented in Gupta et al. (2011), hereafter G11, to check if our age estimator is more sensitive to younger stellar populations. Secondly, we can contrast the local and global properties of hosts to check if there is more information contained in the local environment rather than the more easily studied global properties.

The final analysis of G11 included 206 SNIa and hosts. We looked at the 76 objects that passed our redshift and other quality cuts to use in our method validation. A sample of the data set used for this validation is visible in Table 2.

For each of the 103 hosts where we had local environment photometry, we gathered the SDSS-DR12 model magnitudes for the estimate of the global properties. A sample of this data set can be seen in Table 3.

2.3 Photometry of Nearby Galaxies

The SDSS model magnitudes are unreliable for galaxies with a large angular extent. For the nearby galaxies with distances calibrated with Cepheid variable stars, we used aperture photometry to obtain both the local and global magnitudes. Images of the large galaxies were downloaded from the SDSS DR12 and individual apertures were designed to capture 90% of the combined light in all the filters. After masking out stars projected on the galaxy, the aperture was then applied to the image of each filter. The magnitude was then calibrated using nearby stars. The photometry for these galaxies can be seen in Table 4.

2.4 Supernova Properties

We use the C13 supernova sample to provide light curve properties and Hubble residual information. We use the Malmquist bias corrected distances derived from the best fit cosmology (\text{H}_{0}$$~{}=73.8~{}\text{km s}^{-1}~{}\text{Mpc}^{-1}, $\Omega_{M}=0.24$, and $\Omega_{\Lambda}=0.76$). When needed, we use these values for our assumed cosmology. For the maximum redshift in our sample, the Malmquist bias correction is $\sim 0.01~{}\text{mag}$. However, C13 noted that the stretch correction coefficient, $\alpha$ they found for their full sample was significantly larger than typical and larger than the $\alpha$ derived from their spectroscopically classified sub-sample. After our cuts, we found a significant correlation between the supernova stretch parameters, $x_{1}$, and the C13 Hubble residuals which would likely result in spurious correlations with our host galaxy analysis. We corrected the C13 Hubble residuals using their spectroscopically derived $\alpha$ value of 0.16 and no longer detected a significant correlation between $x_{1}$ and our sample’s Hubble residuals. The resulting Hubble residuals can also be found in Table 1.

For the nearby SNIa used in the Cepheid calibration of SNIa peak luminosity, we obtained light curves from the SNANA database111http://snana.uchicago.edu and fit them using SALT2.4 implemented from the sncosmo222https://sncosmo.readthedocs.io/. The model also corrected for Milky Way dust extinction from the dust maps of Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) via sfdmap333https://github.com/kbarbary/sfdmap.

3 Stellar Population Model

A direct estimate of the age of the stellar population requires a robust model for the observed population. Flexible Stellar Population Synthesis (FSPS) (Conroy et al., 2009; Conroy & Gunn, 2010) takes a star formation history and outputs either a spectrum or a redshifted spectral energy distribution (SED) of the resulting stellar population. The version of FSPS we used (commit ae31b2f from November 2016) uses the MIST isochrones (Dotter, 2016; Choi et al., 2016) and the MILES spectral libraries (Falcón-Barroso et al., 2011).

3.1 FSPS settings

Many of the FSPS parameters were set at their default values, but a number of key settings were adjusted to produce the desired model space.

To control the metallicity of our stellar population, we set zcontinuous = 2. This setting convolves the SSPs (simple stellar populations) with a metallicity distribution function. The metallicity distribution is defined as

[TABLE]

with

[TABLE]

where $z$ is the metallicity in linear units and $z_{\sun}=0.019$. This metallicity distribution is governed by two more FSPS parameters: pmetals and logzsol (i.e. $\log(z/z_{\sun})$). We left pmetals at its default value of $2$ and during the fitting process logzsol was allowed to vary but was marginalized over when the age probability distribution was determined.

The next set of parameters govern the treatment of dust. We used the default power law dust model as explained in Conroy et al. (2009), based off of Charlot & Fall (2000). The attenuation curve of a star, as a function of stellar age, is defined as

[TABLE]

where $\tau_{1}$ and $\tau_{2}$ are the attenuation around a young stellar population and in the ISM respectively. See Charlot & Fall (2000) Figure 1 for a visual representation. In FSPS these two parameters are controlled via the dust1 and dust2 variables. For this analysis, dust1 was set to two times dust2, and dust2 was allowed to vary to match the observations. Conroy et al. (2009) claims good values of dust1 and dust2 are 1.0 and 0.3 respectively, while Charlot & Fall (2000) prefers values of dust1 and dust2 of 1.0 and 0.5 respectively. The allowed range for dust2 in this analysis is explained in Section 4.2 and is consistent with these recommendations.

A few of the host SEDs show an unusual, e.g. SN4019. Adding nebular emission (setting add_neb_emission = True and cloudy_dust = True) adjusts the $r$-band magnitude for a young stellar population and allows the model to match this observed feature. This characteristic is shown to be achievable in the self-consistency validation test number 3, as explained in Section 5.1.

Finally, FSPS outputs the luminosity of $1~{}\text{M}_{\sun}$, so an extra constant, $\delta$, is used to scale the output SED of FSPS to match the observed SEDs.

3.2 Star Formation History

FSPS has many inputs for describing the star formation history of a galaxy. The sfh parameter allows the user to select the functional form of star formation history. G11 used the simple $\tau$-model: the star formation rate is proportional to $e^{-t/\tau}$, with $t$ being the time since the start of star formation and $\tau$ is a free parameter. G11 fit both $\tau$ and the length of star formation history. This is the simplest model, which is important when fitting a small number of free parameters. However, such a simple prescription makes it difficult to create a model with both old stars and recent star formation.

Simha et al. (2014) investigated the ability of several star formation history models to match simulated galaxies. This research looked at the simple $\tau$-model, a linear-exponent model444This is the same as FSPS’s delayed $\tau$-model, sfh = 4., and a four parameter $\tau$-model. A visual comparison is presented in Figures 3–5 of Simha et al.. According to the calculations in Simha et al., the simple $\tau$-model can overestimate the age by $\sim 1\text{--}2~{}\text{Gyr}$, particularly for younger populations. Since we expect some supernovae to explode in young ($\lesssim 2~{}\text{Gyr}$) stellar populations, we decided to use a four parameter $\tau$-model.

The main feature of the four-parameter $\tau$-model (4p$\tau$-model) is that it separates the properties of early and late time star formation. This model can describe a wide range of star formation histories: an early burst, a history dominated by recent star formation, or both an early burst and recent star formation. The 4p$\tau$-model is a piecewise combination of a linear-exponent star formation history then a linearly rising or falling star formation. This model is used by FSPS when sfh = 5. Mathematically the 4p$\tau$-model can be written as

[TABLE]

where $t_{0}$ is when the star formation started, $\tau$ is the e-folding parameter, $t_{i}$ is the star formation transition time, $m_{\text{sf}}$ is the slope of the late time star formation, $\mathscr{R}$ is the ramp function, and $\mathscr{A}(z)$ is the redshift dependant time when the observed light was emitted. Note that the equation above allows negative values of $\Psi$, which is nonphysical. So we add an extra constraint that forces $\Psi(t)$ to be [math] if calculated to be negative. The four variables in the equation ($t_{0}$, $\tau$, $t_{i}$, $m_{\text{sf}}$) are the free parameters that give this model its name. These correspond to the FSPS parameters sf_start, tau, sf_trunc, and sf_slope respectively. Also, this function takes $t$ as the time from the start of the universe. A sample of various star formation histories calculated from the tau model can be seen in Figure 1.

3.3 Calculating Ages

For any set of star formation parameters we determine the average population age. The mass weighted average age is:

[TABLE]

where all of the variables are the same as described in Equation 4 so $t-t_{\text{0}}$ is simply the length of star formation. In the integral, a variable substitution of $t-t_{0}\rightarrow t$ transforms the time zero point from the beginning of the universe to the start of star formation. We would also need to transform $\Psi_{\star}(t)\rightarrow\Psi^{\prime}_{\star}(t)$ so that $\Psi^{\prime}_{\star}$ assumes $t=0$ is the start of star formation. This makes the equation for the mass weighted average age to be

[TABLE]

or the equation used in G11. In this paper, the age of a stellar population will refer to the mass weighted average described here.

The model using the star formation parameters $t_{0}=8.0$, $\tau=0.1$, $t_{i}=12$, and $m_{\text{sf}}=20$, produces a population with an average age of $437~{}\text{Myr}$. This demonstrates that our SFH prescription can generate a dominant young population when SNIa are expected to start exploding. A small $\tau$ is needed to keep the number of old stars from building up over the cosmic time and dominating over the recent linear star formation period. An example of this can be seen by the “old & young burst” star formation history in Figure 1.

4 Determining the most Probable SFH

Using Bayesian statistics and a Markov chain Monte Carlo (MCMC) sampling method we determine the probability of each free parameter described in Section 3: $\log(z/z_{\sun})$, $\tau_{2}$, $\tau$, $t_{0}$, $t_{i}$, $m_{\text{sf}}$, and $\delta$. Data modeling and parameter estimation is often done with Bayesian statistics because it calculates the probability of the model parameters given the observed data by using Bayes’ Theorem:

[TABLE]

where $\theta$ is a given model’s parameters and $D$ represents the observed data. Each probability is:

$P(\theta|D)$

The posterior, which is the probability of the model parameters given the data.

$P(D|\theta)$

The standard likelihood function, $\mathscr{L}(D|\theta)$.

$P(\theta)$

The model prior, which describes what we know about the model before considering the data $D$, such as model parameter limits.

$P(D)$

The model evidence, which in practice amounts to a normalization term.

For MCMC sampling, only relative probabilities are needed, so $P(D)$ is generally ignored and Equation 6 becomes a proportionality, not an equality. It is common to use flat priors, $P(\theta)\propto 1$. In this case Bayes’ Theorem simplifies to a standard frequentist likelihood estimation, $P(\theta|D)\propto\mathscr{L}(D|F)$. More generally, some prior information is used and Bayes’ Theorem becomes

[TABLE]

To find the maximum posterior probability of the model parameters, we need only to know the priors, $P(\theta)$, and the likelihood of the data, $\mathscr{L}(D|\theta)$. A more complete description of Bayesian statistics and MCMC sampling is available in VanderPlas (2014) and Hogg & Foreman-Mackey (2017).

4.1 Likelihood

This method uses a standard log-likelihood function for data with Gaussian uncertainties. Summing over each filter, we compare the observed apparent magnitude ($m_{i}$) with the resulting apparent magnitude from FSPS ($m_{\text{FSPS},i}$) plus a scaling factor ($\delta$) to account for FSPS’s $1~{}\text{M}_{\sun}$ output. Mathematically this is written as:

[TABLE]

with $\sigma_{i}$ as the uncertainty in each $m_{i}$ measurement.

4.2 Priors

For five of the variables we use bounded flat tops as described below:

[TABLE]

Since a flat distribution in a slope parameter preferentially searches the high values space,555A mathematical description is available in VanderPlas (2014, Appendix A); he also has a nice graphical example on his website. the MCMC was performed over $\phi$, the angle the ramp function makes with respect to the $x$-axis; therefore, $m_{\text{sf}}=\arctan(\phi)$. The prior bound above, $-1.520838<\phi<1.520838$, corresponds to $-20<m_{\text{sf}}<20$.

In addition, $\log(z/z_{\sun})$ uses a Gaussian distribution with $\mu=-0.5~{}\text{dex}$ and $\sigma=0.5~{}\text{dex}$ limited to the range of $-2.5<\log(z/z_{\sun})<0.5$. This is a common assumption as seen in Belczynski et al. (2016). Our model reaches a lower metallicity than the grid search used in G11.

For the ISM dust parameter, $\tau_{2}$, we assume a Gaussian prior on the top bounds. The Gaussian distribution is defined by $\mu=0.3$ and $\sigma=0.2$, but only values between [math] and $0.9$ are accepted. This allows for some variability but keeps the values near the $0.3$ and $0.5$ as recommended by Conroy et al. (2009) and Charlot & Fall (2000) respectively.

5 Validation

Following the statistical method described in Section 4, we derive probability distributions for the model parameters defined in Section 3. Using Equation 5 at each step in the MCMC chain, we build a probability distribution function for the age marginalized over metallicity and host galaxy dust.

5.1 Self-consistency

The first validation of this newly developed age estimator was to verify that it was self-consistent, i.e. it could correctly estimate the star formation parameters from an SED generated by FSPS.

Eight different models were used in this test. The initial metallicity, dust, and star formation parameters can be seen in Table 5. These values were put into FSPS, at a redshift of $z=0.05$, to generate observed SEDs. The resulting SEDs, Table 6, were then analyzed with our Bayesian estimator. This set of models produced old populations ($\sim 10.5~{}\text{Gyr}$) and very young populations ($\sim 0.5~{}\text{Gyr}$). They explored the effect of metallicity (Model 5) and dust (Model 6). Model 2 also looked at a “mixed” population with an old burst of star formation and a strong increase of star formation to the present epoch, a stellar population that is not possible with a simpler star formation history.

This method can model a young stellar population, (like Model 3), but is unable to recover a star-burst or mixed populations (Models 8 and 2 respectively) based on SED fitting. The “old & young burst” of Model 2 is not particularly blue and as such, we identify an age of $4.4~{}\text{Gyr}$. We do better with Model 3 where we estimate the correct age of 1.4 Gyr with a small uncertainty.

Figure 2 shows the FSPS output SEDs of the initial input star formation parameters and the best-fit parameters from our analysis. The SEDs are fit very well across all the models used in this self-consistency test. In addition to fitting the SEDs, our analysis was able to approximate the underlying parameters. An example corner plot of the posterior probabilities is displayed in Figure 3, and the the full set of figures is available online.

Fig. Set3. Self-consistency posterior results.

5.2 Comparison to Previous Work

The final validation was to recalculate the global host ages originally presented in G11. A direct comparison between ages from G11 and the results from our analysis can be seen in Figure 4. Most points fit along the one-to-one line with a scatter of around $2~{}\text{Gyr}$ implying that our analysis is consistent with the previous work. However, there are six hosts that our method estimates to be $\lesssim 2~{}\text{Gyr}$, whereas none of the ages in G11 were that young. At $\lesssim 4~{}\text{Gyr}$, our method systematically estimates a younger age. This is because our star formation history model allows for more recent star formation than the simple $\tau$ models permits. In their discussion on this topic, Simha et al. (2014) claimed that the star formation model used in G11 can overestimate young populations by $\sim 2~{}\text{Gyr}$. This overestimation can be seen in Figure 4.

6 Results

Our analysis generates a probability distribution for the age of a stellar population. Probability distributions can be summarized by a median and 68% confidence intervals, especially if the distribution is Gaussian. For non-normal distributions, particularly ones with long tails or multiple modes, the distribution can still be summarized by a median age, but its interpretation is not as clear. To accurately represent the estimated age posterior probabilities shown in Figures 5, 6 and 8, we plotted the results of 100 random samples from both the age and Hubble residual distributions for each SNIa.

We test for correlations between parameters with the Spearman’s rank-order correlation coefficient. This is the non-parametric version of the more common Pearson’s correlation coefficient. There are several differences between Pearson’s and Spearman’s correlations. The most important difference for our study is that the Spearman’s correlation has a high absolute value for any monotonic relationship rather than just linear relationships. This means a linear, exponential, or a single step function would all rank highly with the Spearman’s correlation, but not necessarily with the Pearson’s correlation. Since several previous host galaxy studies have seen steps in Hubble residuals (or more generally sigmoid functions) it is reasonable to use a statistical measurement that is sensitive to these non-linear correlations. See Wall & Jenkins (2012, section 4.2.3) for more information on the Spearman’s correlation.

A large absolute rank-order correlation refers to a tighter association of points, indicated by a small scatter around the correlation. Significance of a Spearman’s correlation can be described by a standard $p$-value, or the probability under the null hypothesis of obtaining a result equal to or more extreme than what was observed. Given a sample size and Spearman’s correlation a $p$-value can be calcuated. If we let 3$\sigma$ be our significance limit, then for our main sample of 103 objects, the Spearman’s critical correlation value would be $\pm 0.30$. That is, there is a 0.2% chance of seeing a Spearman’s correlation value of $>0.3$ or $<-0.3$ from our data set assuming no underlying correlation.

6.1 Global Environments

First, we compare the Hubble residuals with ages derived from the global photometry of the hosts using C13 (Table 3) sample. The comparison is presented in Figure 5. This data set has a low 2.1$\sigma$ correlation, as defined by the Spearman’s value of -0.23. In addition, the distribution in Hubble residual-age space appears to show a distinctive “step” between 7 and $8~{}\text{Gyr}{}$.

6.2 Local Environments

Finally, we compare the Hubble residuals versus average local environment age for the data set derived from C13 (Table 1). These results are presented in Figure 6. Using an age derived from the local environment only slightly changes the Spearman’s correlation between these two parameters, But this correlation, -0.21, only has a 1.8$\sigma$ significance. The overall age distribution and apparent “step” at $\sim 8~{}\text{Gyr}{}$ are not significantly changed by switching to a local environment analysis. At first glance, the local age does not appear to contain any additional information not already present in the global age. Jones et al. (2018) found a marginally significant difference between their global and local analyses.

7 Analysis and Discussion

7.1 Comparison Between Local and Global Ages

The global age estimated for spiral galaxies is an average of several stellar populations. The prompt component of SNIa are expected to be found in the youngest regions of a galaxy. So we may see that the population age at the supernova location would tend to be younger than the global age of the galaxy. In Figure 7 we show the difference between the estimated global and local ages. For massive galaxies with local ages less than 4 Gyr there is a tendency for the supernova to be located in a younger than average spot in the host. The effect is less apparent for low mass galaxies, but that is likely because the size of the region measured for the local environment is a large fraction of the size of a small galaxy.

Between local ages of 4 to 8 Gyr the difference between global and local age estimates show a large scatter with no apparent trend. Beyond 8 Gyr the scatter between the local and global ages is reduced, probably because the stellar population in ellipticals is fairly uniform. In these older hosts the global age estimates tend to be 1 to 2 Gyr younger than the local ages. The reason for this difference is not clear, but it may be due to activity at the center of ellipticals contaminating the stellar colors.

Although we see no statistical difference between the use of local and global ages when comparing with Hubble residuals, our results suggest that using local photometry to characterize the supernova environment has some benefit over the global average. For example, the SED measured local to the supernova can indicate a population that is a factor of two younger than the global average in large star-forming galaxies. When feasible, the measurement of the local environment, particularly of younger populations, provides a more accurate representation of the progenitor age than simple averaging the light from the host.

7.2 Investigating the Age Step

The Hubble residuals in the C13 sample has a mono-tomically decreasing trend (at 2.1$\sigma$), that appears to be a break or step at an age of $\sim 8~{}\text{Gyr}$. Figure 8 plots the same data as Figure 6 but this time splits the data into two age bins: $\leq 8~{}\text{Gyr}$ and $>8~{}\text{Gyr}$. Both age ranges show very small Spearman’s correlations ($-0.07$ and $-0.09$ respectively) that are consistent with a flat distribution. The correlation decreasing in significance when the data set is split in two, implies that the monotonic function seen in Figure 6 is really a step-like function with a transition at $\sim 8~{}\text{Gyr}$. This function may have a transition width making it more like a continuous sigmoid function with a transition faster than our age resolution. In both the local and global age measurements, we find a significant age step in Hubble residuals.

The amplitude of the jump between these two populations split at $\sim$ $8~{}\text{Gyr}$ appears to be more than 0.1 mag. The younger population has a mean Hubble residual just over zero ($0.005\pm 0.018~{}\text{mag}$), whereas the older population has a mean Hubble residual of $-0.109\pm 0.035~{}\text{mag}$. This makes the resulting difference in the means equal to $0.114\pm 0.039~{}\text{mag}$ ($2.9\sigma$). This is almost two times larger than the commonly used mass step of $0.06~{}\text{mag}$. A step of this size may affect precision cosmological measurements since the fraction of each subpopulation is expected to change with redshift. For example, one would not expect to find stellar populations as old a $8~{}\text{Gyr}$ at redshifts greater than $\sim 1$ given the current standard cosmology. Such an age step may also impact local measurements of the Hubble constant as the peak luminosity of SNIa have so far been calibrated with Cepheid variables found exclusively in star forming galaxies. With the age step being detected at a $3\sigma$ significance, Figure 8 suggests a need for an additional SNIa luminosity correction based on their host’s local or global stellar age.

7.3 Age as the Cause of the Mass Step

Childress et al. (2014) argues for a link between host stellar mass and the delay time between stellar formation and SNIa explosion, with results summarized in their Figure 4 where the SNIa progenitor age distribution is divided into host galaxy stellar mass bins. They find, with reasonable assumptions of star formation histories and SNIa delay times, that there is a natural division in host mass and age between prompt and “tardy” SNIa. Prompt SNIa occur in lower mass galaxies with ages $\lesssim 6~{}\text{Gyr}$ while tardy events continue in high mass galaxies with ages $\gtrsim 6~{}\text{Gyr}$. Projecting their model onto the host mass axis results in an overlapping distribution of prompt/tardy explosions with a transition near 10${}^{10.5}~{}\text{M}_{\odot}$. This transition in progenitor age may correspond to the mass step observed at 10${}^{10}~{}\text{M}_{\odot}$. Projection of the Childress et al. (2014) model onto the age axis provides a clean separation between the prompt and tardy SNIa with a dearth of events between 4 and 8 Gyr. Our age estimates are not consistent with a deficit of supernovae exploding in that range, but we do see a shift in the SNIa light curve properties around a stellar age of about $8~{}\text{Gyr}$.

We have estimated the host galaxy masses in our sample to test if our measurements agree with the prediction of Figure 4 in Childress et al. (2014). We use the kcorrect code (v4_3)666Available through http://kcorrect.org or http://github.com/blanton144/kcorrect. described in Blanton & Roweis (2007). This code utilizes spectral fitting templates based on the stellar population synthesis models of Bruzual & Charlot (2003), which are calculated using the Chabrier (2003) initial mass function. We input the SDSS model magnitudes and the pipeline redshifts for each galaxy (presented in Table 3) to calculate the $k$-corrections and stellar mass-to-light ratios. The stellar masses are output in units of M${}_{\odot}h^{-1}$, and we convert them to units of M*⊙* using the C13 cosmology described in Section 2.4. The uncertainties on the stellar masses are approximately $\pm$0.3 dex. The resulting estimated stellar masses are reported in Table 7.

Figure 9 shows the distribution of the local SNIa age versus stellar mass for our sample. For our sample, the distribution of hosts in age-stellar mass space shows similarities and differences with that predicted in Childress et al. (2014). The observations do show young hosts extending to low stellar masses reproducing the backward “L” seen in the Childress et al. (2014) simulation. This indication of “cosmic downsizing” is not as pronounced in our data as we are probably still overestimating the ages of the extremely young populations. Most notably missing, is that the predicted bi-modal age feature is not present. An island of young galaxies is expected between $0.5-1.0~{}\text{Gyr}$ and that is not seen in Figure 9. Our age distribution is relatively flat from 2 Gyr to 5 Gyr and then declines out to 11 Gyr. This does not match the predicted distribution in Figure 4 of Childress et al. (2014), and there is no clear peak of old hosts. The fact that our sample of supernovae extends to a redshift of 0.2 may contribute to the lack of a clear peak in old age. Childress et al. does predict that the old-age peak decreases in height and age as the data is taken at higher redshifts simply due to the finite age of the universe.

7.4 Principal Component Analysis

Light curve shape and color of SNIa have been shown to be strongly correlated with peak luminosity. Other modest trends with population age, host mass, and gas metallicity have been reported and here we have identified a rather significant jump in Hubble residuals with local and global stellar age. Because these environmental observables are highly correlated with each other, it is interesting to ask if there is a linear combination of observables that have a clear and significant correlation with SNIa Hubble residual. As an initial test, we performed a principal component analysis (PCA) on our data set.

PCA is a linear algebra tool that transforms the basis of a matrix of data to orthogonal axes where the new axes are maximally aligned and sorted with the intrinsic scatter of the data.777For more information, see Ivezić et al. (2014, section 7.3), Wall & Jenkins (2012, section 4.5), or visit https://towardsdatascience.com/pca-using-python-scikit-learn-e653f8989e60 Therefore the first principal component contains the most amount of information and the last principal component contains the least. Thus, it is possible to reduce the dimensions of a problem by retaining only the principal components that comprise most of the variance (or relative information) in the original data set. In addition, the first principal component will identify a linear combination of the original variables that accounts for most of the variance in the data. Interpretation of PCA is difficult because the results are sensitive to noise, specific parameterizations, normalization, and other implementation details. However, PCA can still be useful to understand how multiple parameters work together within a data set. In the search for a trend in Hubble residuals versus SNIa parameters, we applied PCA on the parameters of SALT2 stretch ($x_{1}$) and color ($c$), as well as host mass and age. We calculated PCA coefficients first using local ages and then with global ages. Hubble residuals were not included as a PCA parameter. Only after the PCA did we search for correlations between the resulting principal components and the Hubble residuals.

Before running PCA, we normalized all parameters by removing the mean and scaling to unit variance resulting in normalized input parameters ($x_{1}^{\prime}$, $c^{\prime}$, $m^{\prime}$, and $a^{\prime}$). An example, for $x_{1}$, is defined as:

[TABLE]

The means and standard deviations used in the normalization process can be seen in Table 8.

From these four observables, PCA yields four principal components (PCi). Table 9 shows the linear combination of the observable variables that make up each principal component as well as the explained variance. Substituting global age for local age resulted in some minor differences between the PCA coefficients and variance, but the overall conclusions are very similar between the two analyses. PC2 accounts for a quarter of the variance and is dominated by the SALT2 color parameter. The other two PCA components contribute only a small portion of the variance.

As a first analysis, it is important to see if there are any correlations between Hubble residual and these SNIa-host galaxy prinicpal components. Looking at the PCA done with local age, PC3 and PC4 versus Hubble residual have extremely low Spearman’s coefficients of 0.21 and $-0.15$ respectively. PC2, dominated by the SALT2 color term, when compared with Hubble residual is also a scatter plot with a low Spearman’s correlation of $-0.06$. Surprisingly there is an increase in the scatter at the high PC2 (high color) domain. Using the PCA done with the global age, PC2, PC3 and PC4 have similar Spearman’s correlation coefficients. Figures showing the relationship between Hubble residual and each principal component (4 per analysis method, 8 total figures) are available as a Figure Set in the online journal.

Fig. Set10. Hubble residuals versus principal component

Figure 10 shows the relationship between Hubble residuals and PC1 using the estimated local age. The correlation when PCA is applied using the global host age can be seen in Figure 11 (a). PC1 with either local or global ages shows a very strong correlations with Hubble residuals. The first principle component, using local age estimates, is defined as:

[TABLE]

or it can be approximated by ignoring the color term since it barely contributes to PC1.

In PC1, the SALT2 color coefficient is very small, implying supernova color is not a strong contributor to this Hubble residual correlation. The mass, age, and stretch parameters in our PCA analysis are similar in amplitude and likely are of similar importance in any further standardization of SNIa distances. The SALT2 stretch parameter, $x_{1}$, has a surprising large contribution to PC1 given that the SNIa have already been corrected for light curve shape. Figure 10 suggests a correlation between stretch and Hubble residual, but this correlation was removed at the start of the analysis. Instead, PC1 shows that the value of $x_{1}$ is related to mass and age. Attempting correct for the stretch-luminosity relationship requires including host properties as all three have an influence on Hubble residuals. Our results suggest that the $\alpha$ parameter derived from the SALT2 fit is not ideal, because the affects of host mass and population age are not distributed uniformly with stretch.

The Hubble residual-PC1 correlation in Figure 10 has a Spearman’s correlation of 0.44. This trend is highly significant with a p-value of $3.86\times 10^{-6}$ corresponding to a $4.7\sigma$ significance. The trend in Hubble residual versus age has a significance of only 2.1$\sigma$ (Figure 5), while including mass and stretch greatly increases the significance. This strong correlation suggests that host properties influence the luminosities of SNIa beyond the currently applied corrections.

When using global ages, the underling trend has a Spearman’s correlation of 0.40 or a $4.0\sigma$ significance. The Spearman’s correlation with PC1 from local ages is slightly more significant than for global ages suggesting that measurements near the event provide some improvement in correlating supernova with environment.

The best-fit linear regression of this trend is

[TABLE]

with uncertainties in the slope and intercept as $\pm 0.011~{}\text{mag}$ and $\pm 0.016~{}\text{mag}$ respectively. This trend reduces the $1\sigma$ scatter in Hubble residual from $0.17~{}\text{mag}$ to $0.15~{}\text{mag}$. The best-fit linear regression using global ages values is ${\rm HR}=0.045\times{\rm PC}_{1}-0.012$ with uncertainties in the slope and intercept as $\pm 0.011~{}\text{mag}$ and $\pm 0.016~{}\text{mag}$ respectively. This reduces the $1\sigma$ scatter in the Hubble residual to $0.16~{}\text{mag}$. As a note, the PCA does not maximize this correlation or minimize the $\chi^{2}$ parameter in a linear regression. For example, reversing the sign of the $x_{1}$ coefficient nearly doubles the slope of the correlation with Hubble residual as seen in Figure 11 (b) and significantly reduces the $\chi^{2}$ parameter of a linear fit. Understanding the best use of these coefficients will be part of future work already in preparation.

Since the age and mass coefficients make up nearly two-thirds of the correlation amplitude, we see that a change of $\pm 2\sigma$ in the normalized mass and age values results in a 0.24 mag shift in SNIa brightness. Cosmic downsizing suggests these two parameters are correlated since young galaxies tend to be small while typical old hosts are massive in the current epoch. Low-mass, young galaxies tend to be metal deficient while old, massive galaxies can be metal rich. Thus, this combination of age and mass may indicate progenitor metallicity influences the peak luminosities of SNIa.

The correlation between Hubble residual and PC1 is very significant, and it is unlikely that PCA could generate such a strong correlation from a random distribution. To test the probability that this correlation is caused by chance, we applied a bootstrap style method for hypothesis testing. We took our PCA input matrix ($x_{1}$, $c$, host mass, local environment age) and shuffled the order along each parameter, creating 103 “new” SNIa. There was no cross shuffling, so a stretch always stays a stretch, it just corresponds to a different SNIa. We applied PCA on each shuffled sample and tested for any correlations with Hubble residual. After 32,000 runs, a $4\sigma$ test, the maximum Spearman’s correlation coefficient between Hubble residual and any principal component was 0.43 while the measured Spearman’s test of the non-shuffled values was 0.44. The distribution of the bootstrapped Spearman’s correlation coefficients is shown in Figure 13. This bootstrap style analysis shows that the false positive Spearman’s correlation follows a Gaussian distribution with a $\sigma\approx 0.1$. This is the expected distribution for Spearman’s correlation coefficients for a data set of $N\approx 100$. It is exceedingly unlikely that the correlation between Hubble residual and PC1 found in Figure 10 would appear at random. As a result, this luminosity-stretch-mass-age relationship is the most significant systematic seen between calibrated SNIa and host-galaxy environment.

7.5 Correcting for the Hubble Residual-PC1 Correlation

Since our PC1 strongly correlates with Hubble residuals, it is reasonable to consider a modification to the Tripp formula that is used to correct SNIa peak luminosities. The new correction coefficients would be the multiplication of the PCA coefficients and the slope of the Hubble residual-PC1 trend. Except for the SALT color, the individual components making up PC1 have significant weights, so are included in this modified equation. The stretch parameter in PC1 can be grouped with the original Tripp coefficient, leaving a new term with just host properties. Since PC1 has a positive correlation with Hubble residual, the coefficient should come in with a negative sign. From these results, and following the example of others (e.g. Moreno-Raya et al., 2016b), we propose a change in the distance modulus corrections performed by SALT2 by modifying the equation to include PC1. This new equation would be

[TABLE]

where $\alpha^{\prime}\approx\gamma\approx 0.03~{}\text{mag}$. This uses both the approximate form of PC1 and encompasses the very slight renormalization of $x_{1}$ into $x_{1}^{\prime}$. More research is needed to accurately determine $\alpha^{\prime}$ and $\gamma$.

This modified distance correction formula now includes a host galaxy stellar mass correction and a host age term. Childress et al. (2014) attempted to explain the mass step as an age dependency and we see that both appear to have an impact on SNIa luminosities. Our result is bolstered by Jones et al. (2018) who showed that stellar population color effects are still present after stellar mass corrections. Stellar mass plus population age is similar to the Mannucci relationship (Mannucci et al., 2010) that connects galaxy mass and star-formation rate with gas phase oxygen abundance. Our combined age plus mass parameters work in the same sense as the Mannucci relation and maybe a proxy for metallicity. Both Hayden et al. (2013) and Moreno-Raya et al. (2016b) have indicated that metallicty is the underlying galactic variable influencing Hubble residuals and our results support this suggestion.

Kasen et al. (2009) shows that varying the metallicity of SNIa progenitors will shift the slope of the Phillips relationship. To investigate this we assume our SNIa have a ${\rm M_{peak}}=-19$. We then add in the Phillips relationship (with $\alpha=0.16$) and our Hubble residual-PC1 relationship. The result is seen in Figure 14. Since, fast decliners are preferentially in larger and older hosts, you see the host galaxy effect on the Phillips relationship is stretch dependant. Put another way, the slope of the Phillips relationship, $\alpha$, is dependant on host galaxy properties of the sample. These results show that these parameters are interrelated and that there is likely a multi-dimensional relationship between peak absolute magnitude, decline rate, and progenitor metallicity. Research into this multi-dimensional relationship will be a part of a future work already in preparation.

8 The Affect on $\text{H}_{0}$

The most precise SNIa measurement of the local value of $\text{H}_{0}$ was presented in Riess et al. (2016, 2018). Riess et al. rebuilt the distance ladder connecting geometric distances, Cepheid distances and finally SNIa in the Hubble flow. These works take into account the host galaxy mass step, but this correction has a minor impact (0.7%, Riess et al., 2016) on the resulting $\text{H}_{0}$. Our trend in PC1 includes stellar age that has the potential for a significant impact on $\text{H}_{0}$ since Cepheid hosts tend to have a significantly younger stellar population than the average Hubble flow galaxy.

To test if the correlation in the PC1 parameter could affect $\text{H}_{0}$ , we estimated the local age, global age and stellar mass of the SNIa hosts that have distances calibrated using Cepheid variables. We then compare the average PC1 parameter found for the Cepheid sample with the average for a Hubble flow sample. As a representative Hubble flow data set, we use our analysis of the C13 galaxies. The difference in the average PC1 values for the two sets of hosts, multiplied by the slope in Equation 12 provides an estimate of the shift in peak absolute magnitude, $\Delta M$ between the calibration sample and the Hubble flow. The fractional error in distance, and therefore $\text{H}_{0}$, due to the differences in age, mass, and stretch between the two samples is then

[TABLE]

SDSS imaging was available for 14 of the 19 SNIa hosts presented in Table 1 of Riess et al. (2016). The others were outside of the SDSS footprint and without SDSS-$u$ photometry, and therefore we were not able to apply our age estimator in a consistent way. For this calibration sample, we followed the same procedures to get the local environment age, global age, and host stellar mass. These values, along with the full list of SALT2 parameters can be found in Table 11. A visual comparison of these parameter’s distributions can be seen in Figure 15.

The individual parameters of PC1 have only mild differences between the SNIa calibration sample and the Hubble flow hosts. The most extreme fast decliners ($x_{1}<-2$) are only in the Hubble flow sample because they are preferentially found in passive galaxies. The mass distributions are significantly different, with a low-mass tail in the calibration sample. The Hubble flow sample has a wider distribution of local ages than the calibration set, same with the global age. What is interesting about PC1 is that all these slight differences work in the same direction. The lower stretch, higher mass and older ages in the Hubble flow combine for a significantly lower average PC1 when compared with the larger stretch, lower mass, younger hosts from the calibration sample.

For the calibration sample, we calculated a PC1 from both the local and global age PCA methodology, PC1,local and PC1,global respectively. These values can also be seen in Table 11. A comparison of these two PC1 distributions are shown in Figures 16 and 17.

Using the local age, the calibration sample in Figure 16 has a mean PC1,local of $0.41\pm 0.15$ ($2.7\sigma$), while the Hubble flow PC1 mean is defined as zero. To account for the sample sizes, we performed the Kolmogorov-Smirnov test (KS-test) resulting in a 5.5% chance ($1.9\sigma$) that these samples are drawn from a common distribution and a Mann-Whitney U test indicates that the calibration sample has a higher mean at a $1.2\sigma$ level. Using the best fit of the correlation between Hubble residual and PC1, this corresponds to a difference of $0.021~{}\text{mag}$.

Rerunning this analysis with the global age PCA normalization and methodology, there is a $0.63\pm 0.14$ ($4.6\sigma$) shift in the PC1,global means. The resulting KS-test says there is a 0.8% chance ($2.7\sigma$) that these two samples are drawn from a common distribution. A Mann-Whitney U test indicates that the calibration sample has a higher mean at a $2.3\sigma$ significance. Applying this shift to the trend seen in Figure 11 (a), there is a shift in peak luminosity of $0.028~{}\text{mag}$.

These shifts in peak luminosity produces a 1.0% or a 1.3% effect on $\text{H}_{0}$ respectively. This is about twice the size of the already accounted for mass step, but is less than the current $1\sigma$ uncertainty in $\text{H}_{0}$ (2.3%). A large SNIa systematic effect was found but it had a minimal effect on $\text{H}_{0}$. For this effect to relieve the full $3.8\sigma$ tension (Riess et al., 2018), these two samples would need to have a PC1 shift of $\sim 3.5$, about 6 times larger than currently seen. The SNIa systematic found in this paper is extremely unlikely to fully resolve the $\text{H}_{0}$ tension.

9 Conclusion

Host galaxy properties have an effect on the absolute magnitude of SNIa. Host galaxy stellar mass, age, and metallicity have all been shown to be a secondary correction to the Phillips relation. Using a Bayesian method to estimate the age, we were able to look at how Hubble residuals of SNIa correlate with the mass weighted average age for both the local environment and the galaxy as a whole. This method is better at correctly estimating younger populations than previous methods. A 2.1$\sigma$ significant correlation between Hubble residual and age was seen. This correlation may be an age step of $0.114\pm 0.039~{}\text{mag}$ at $\sim{8~{}\text{Gyr}}$. This step is nearly twice the size of the currently used mass step.

Running this analysis on both the local environment and the galaxy as a whole showed that the local age did not show any stronger of a systematic than the global age, but as expected, the local age is younger than the global age of the host for SNIa in young populations.

We are unable to completely replicate the predictions of Childress et al. (2014) for the distribution of SNIa hosts in the space of age versus stellar mass. The bi-modal SNIa age distribution is not present in the data set derived from C13, but the general trends of old galaxies with high stellar masses and a tail of low-mass young galaxies are seen in the data.

Using PCA on the two SALT2 parameters, host stellar mass, and local environment age, we see a very significant correlation (0.44) between Hubble residual and the first principal component ${\rm PC}_{1}=0.56x_{1}-0.10c-0.54m^{\prime}-0.63a^{\prime}$. This trend was fit with a slope of $0.051\pm 0.011~{}\text{mag}$. The mixture of parameters making up PC1 suggests that to understand the luminosity variations in SNIa and to properly correct for them requires simultaneous knowledge of their host and supernova properties.

As a result of this significant trend, PC1 should be used as part of an an updated light curve fitter.

The dominant components of PC1 are stretch, mass, and age. Using the Mannucci relationship, PC1 may be implying that $\alpha$ has a metallicity dependence. A theoretical case for this was already made by Kasen et al. (2009), and the observational trends found in this our work are able to reproduce this predicted effect.

A correlation of this magnitude could have major effects on the precision measurement of $\text{H}_{0}$. Looking at the difference in the calibration and Hubble flow sample, we see that these data sets have a meaningful difference in PC1. Using the PCA methodology and normalization from the global age analysis, there is a shift in PC1,global of $0.63\pm 0.14$ (4.6$\sigma$). In addition, a KS-test shows that they are drawn from different underling populations at the $2.7\sigma$ significance level, and a Mann-Whitney U test says the calibration sample has a higher mean at a $2.3\sigma$ significance. Similar differences in the mean were seen using the local age. This difference between these two samples would correspond to a shift in SNIa peak absolute magnitude of $0.028~{}\text{mag}$, or a 1.3% shift in $\text{H}_{0}$. With at most a $\sim 0.5\sigma$ effect on $\text{H}_{0}$, this is extremely unlikely to relieve the full $3.8\sigma$ tension between the most recent measurements of the CMB from the Planck collaboration.

A major systematic in SNIa was discovered, but it had only a small effect on $\text{H}_{0}$. This correction should be further investigated and applied to SNIa used in cosmological studies. Moreover, it appears that even a large SNIa systematic cannot fully relieve the tension between the local and CMB measurements of $\text{H}_{0}$.

The authors would like to thank Chris Wotta and Eric Bechter who provided invaluable feedback on improvements to the analysis code’s design and Erika Holmbeck for suggestions on the figures. In addition, thank you to David Rubin, Daniel Foreman-Mackey, David Jones, Adam Riess, and the anonymous referee for their informative comments. Funding for this research was in part by David Eartly and the Lennox Fellowship. This research was supported in part by the Notre Dame Center for Research Computing. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington. SDSS