Spectrophotometric templates for core collapse supernovae and their application in simulations of time-domain surveys

TL;DR

This paper introduces data-driven spectral templates for core collapse supernovae, enabling more accurate simulations of their populations in time-domain sky surveys, which is crucial for survey design and analysis.

Contribution

The paper provides a comprehensive set of spectral templates for various core collapse supernova types, extending into UV, with open-source code for broad application in simulations.

Findings

Templates cover hydrogen-rich and stripped-envelope supernovae.

Templates extend into near-ultraviolet wavelengths.

Demonstrated integration with SNANA for supernova simulations.

Abstract

The design and analysis of time-domain sky surveys requires the ability to simulate accurately realistic populations of core collapse supernova (SN) events. We present a set of spectral time-series templates designed for this purpose, for both hydrogen-rich (type II, IIn, IIb) and stripped envelope (types Ib, Ic, Ic-BL) core collapse supernovae. We use photometric and spectroscopic data for 67 core collapse supernovae from the literature, and for each generate a time-series spectral template. The techniques used to build the templates are fully data-driven with no assumption of any parametric form or model for the light curves. The template-building code is open-source, and can be applied to any transient for which well-sampled multi-band photometry and multiple spectroscopic observations are available. We extend these spectral templates into the near-ultraviolet to $\lambda \lambda \sim 1600 \AA$ using observer-frame ultraviolet photometry. We also provide a set of templates corrected for host galaxy dust extinction, and provide a set of luminosity functions that can be used with our spectral templates in simulations. We give an example of how these templates can be used by integrating them within the popular SN simulation package $\textsc{snana}$, and simulating core collapse supernovae in photometrically-selected cosmological type Ia supernova samples, prone to contamination from core collapse events.

Tables

Table 1. Table 1: Summary of luminosity functions and relative rates for core collapse SNe

SN type	LFs from	Revised LFs from	LFs from	LFs adjusted from	Ratesd
	Li et al. (2011)a	Li et al. (2011) b	Richardson et al. (2014)c	Jones et al. (2017)
II	-	-15.97(1.31)	-	-	64.9
IIL	-17.44(0.64)	-	-17.98(0.86)f	-18.28(0.45)	7.9
IIP	-15.66(1.23)	-	-16.75(0.98)f	-16.67(1.08)	57.0
IIb	-16.65(1.30)	-16.69(1.38)	-16.99(0.92)	-16.69(1.99)	10.9
IIn	-16.86(1.61)	-17.90(0.95)	-18.53(1.36)	-17.66(1.08)	4.7
Ic	-16.04(1.28)	-16.75(0.97)	-	-17.44(0.66)	7.5
Ib	-17.01(0.41)	-16.07(1.34)	-17.45(1.12)	-18.26(0.15)	10.8
Ibc-pec	-15.50(1.21)	-	-	-	-
Ic/Ic-pec/Ic-BL	-	-16.79(0.95)	-17.66(1.18)	-	8.6
Ic-BL	-	-	-	-	1.1

Table 2. Table 2: The core collapse SNe included in this work.

Name	Type	Redshift	Peak MB	Optical data	Near-UV data	Numb. of	Ref.(c)	$E{\left(B-V\right)}_{host}^{\left(a\right)}$	Ref.
						Spectra			Host(c)
SN2013by	II	0.00359	-18.36	$BVugri$	u,v,uvw1,uvm2,uvw2	7	(68),(69),(42)	0	(69)
ASASSN15oz	II	0.00693	-18.19	$UBVRIgri$	u,b,v,uvw1,uvm2,uvw2	13	(2)	0	(2)
SN2014G	II	0.00563	-18.10	$UBVRI$	v,uvw1,uvm2,uvw2	17	(77),(71),(78)	0.24	(78)
SN2007od	II	0.00586	-17.93	$UBVRIri$	u,b,v,uvw1,uvm2,uvw2	13	(29)	0	(29)
SN2013ej	II	0.00219	-17.76	$UBVRIgriz$	u,b,v,uvw1,uvm2,uvw2	28	(73)	0.049	(94)
SN2013fs	II	0.01185	-17.73	$BVRIr$	u,uvw1,uvm2,uvw2	23	(74),(63),(42),(75)	0.015	(75)
SN2008bj	II	0.019	-17.48	$UBVri$		16	(31)	0.081†	(90)
SN2016X	II	0.00441	-17.38	$UBVRIgri$	u,b,v,uvw1,uvm2,uvw2	33	(79)	0.02	(79)
SN2009bw	II	0.00382	-17.14	$UBVRI$	uvm2,uvw2	16	(39)	0.08	(39)
SN2013ab	II	0.00532	-17.14	$BVRIgri$	u,b,uvw1,uvm2,uvw2	22	(64),(65)	0.02	(65)
SN2004et	II	0.00105	-17.09	$UBVRI$		28	(16)	0.07	(87)
SN2012aw	II	0.0026	-17.05	$UBVRIugriz$	u,b,v,uvw1,uvm2,uvw2	30	(61),(62),(63)	0.028	(62)
SN2008in	II	0.00522	-16.96	$BVri$	u,b,v,uvw1,uvm2,uvw2	10	(36),(31)	0.076	(36)
SN2009dd	II	0.00246	-16.81	$UBVRIri$		6	(30),(31)	0.433	(30)
SN1999em	II	0.00239	-16.69	$UBVRI$		27	(12)	0.077	(86)
SN2009kr	II	0.0065	-16.57	$BVRI$	u,b,v,uvw1,uvm2,uvw2	6	(44)	0	(44)
ASASSN14jb	II	0.00603	-16.11	$BVgri$	u,b,uvw1,uvm2,uvw2	10	(1)	0	(1)
SN2012A	II	0.00231	-16.10	$UBVRI$		39	(57)	0.012	(57)
SN2009ib	II	0.00435	-15.86	$BVRIugriz$		9	(40)	0.131	(40)
SN2013am	II	0.00269	-15.80	$UBVRIi$	u,b,uvw1,uvm2,uvw2	14	(66),(67)	0.55	(67)
SN2009N	II	0.00345	-15.03	$BVRgriz$	u,b,v,uvw1,uvm2,uvw2	22	(37)	0.113	(37)
SN2005cs	II	0.00154	-14.90	$UBVRIz$	b,uvw1,uvm2,uvw2	18	(23)	0.015	(89)
SN2016bkv	II	0.002	-14.89	$UBVRIgri$		15	(80),(81)	0	(81)
SN1987A	87A-like	0.00001	-14.33	$UBVRI$		18	(3)	0.17	(85)
SN2011fu	IIb	0.0185	-18.62	$UBVRIz$		28	(52),(53),(21)	0.015	(53)
SN2006T	IIb	0.0081	-17.72	$BVugri$		6	(8),(19),(18)	0.277 (0.048)(b)	(88)
SN2008ax	IIb	0.0019	-16.65	$UBVRIugriz$	u,b,v,uvw1,uvw2	42	(34),(8),(24)	0.28	(91)
SN2016gkg	IIb	0.00492	-16.61	$BVRIgr$	u,b,v,uvw1,uvm2,uvw2	17	(82)	0.09	(95)
SN2008aq	IIb	0.00797	-16.60	$BVugri$	u,b,v,uvw1,uvw2	7	(8),(18)	0 (0)(b)	(88)
SN2011ei	IIb	0.0093	-16.60	$BVRI$	u,b,v,uvw1,uvw2	13	(50),(51)	0.18	(51)
SN2011dh	IIb	0.0015	-15.99	$UBVRIg$		39	(48),(49)	0.05	(49)
SN2013df	IIb	0.00239	-15.90	$UBVRI$	u,b,v,uvw1,uvm2,uvw2	12	(70),(71),(72),(21)	0.08	(72)
SN2011hs	IIb	0.00571	-15.71	$BVRIugriz$	u,b,v,uvw1,uvw2	13	(54)	0.16	(54)
SN2008bo	IIb	0.005	-15.46	$BVri$	u,uvw1,uvw2	13	(8),(24),(21)	0.033	(90)
SN1993J	IIb	-0.0001	-13.22	$UBVRI$		50	(4),(5)	0.1	(5)
SN2008fq	IIn	0.0106	-18.89	$BVugri$		10	(35),(25)	0.46	(25)
SN2010al	IIn	0.017	-18.75	$BVRIuri$	u,b,v,uvw1,uvm2,uvw2	11	(45),(46),(31),(21)	0.053	(46)
SN2007pk	IIn	0.01665	-18.71	$BVri$	u,b,v,uvw1,uvm2,uvw2	9	(30),(31)	0.081†	(91)
SN2009ip	IIn	0.00594	-18.33	$BVRI$	u,uvw1,uvm2,uvw2	55	(41),(42)	0.01	(91)
SN2006aa	IIn	0.0207	-17.61	$BVugri$		6	(25)	0.081†	(90)
SN2011ht	IIn	0.0036	-16.98	$BVRI$	u,b,v,uvw1,uvm2,uvw2	12	(55),(56)	0.061	(56)
SN2005bf	Ib	0.0189	-18.44	$BVugri$		26	(22),(8),(18),(21)	0.1	(22)
SN2007uy	Ib	0.007	-18.34	$BVri$	u,b,v,uvw1,uvm2,uvw2	11	(8),(24),(21)	0.63	(93)
SN2005hg	Ib	0.021	-18.02	$UBVri$		20	(8),(24),(21)	0.2235†	(90)
SN2009iz	Ib	0.0142	-17.69	$BVuri$	uvw1	10	(8),(24),(21)	0.2235†	(90)
SN2009jf	Ib	0.008	-17.68	$UBVRIuugrriiz$		31	(43),(8),(24),(21)	0.03	(91)
SN2004gv	Ib	0.02	-17.39	$BVugri$		5	(8),(18),(21)	0.053 (0.118)(b)	(88)
SN2012au	Ib	0.0045	-17.39	$BVRIz$	u,b,v,uvw1,uvm2,uvw2	5	(59),(60),(21)	0	(60)
SN2006ep	Ib	0.015	-17.16	$BVugri$		7	(8),(18),(21)	0.233 (0.226)(b)	(88)
SN2004gq	Ib	0.00647	-16.97	$BVugri$		24	(8),(19),(18)	0.067 (0.143)(b)	(88)
SN1999dn	Ib	0.0093	-16.62	$UBVRI$		13	(11)	0.048	(11)
iPTF13bvn	Ib	0.00449	-16.54	$UBVRIgriz$	u,b,v,uvw1,uvw2	25	(83),(84)	0.17	(84)
SN2008D	Ib	0.007	-16.47	$BVri$	u,v,uvw1	20	(33),(8),(24),(21)	0.6	(91)
SN2007Y	Ib	0.0046	-16.32	$BVugri$	u,b,v,uvw1,uvm2,uvw2	6	(27),(18),(21)	0.090 (0.318)(b)	(27)

Table 3. Table 3: continued

Name	Type	Redshift	Peak MB	Optical data	Near-UV data	Numb. of	Ref.(c)	$E{\left(B-V\right)}_{host}^{\left(a\right)}$	Ref.
						Spectra			Host(c)
SN2004aw	Ic	0.0159	-17.66	$UBVRI$		23	(15)	0.35	(15)
SN1994I	Ic	0.0015	-17.05	$UBVRI$		27	(6),(7),(8)	0.42	(7)
SN2004fe	Ic	0.018	-17.05	$BVugri$		9	(17),(8),(18)	0.000 (0.052)(b)	(88)
SN2004gt	Ic	0.00548	-16.66	$BVugri$		13	(20),(8),(18),(21)	0.237	(88)
SN2013ge	Ic	0.0055	-16.57	$VRIri$	u,b,v,uvw1,uvm2,uvw2	18	(76)	0.047	(76)
SN2007gr	Ic	0.0017	-16.20	$UBVRIri$		36	(28),(8),(24),(21)	0.03	(92)
SN2011bm	Ic	0.0022	-13.10	$UBVRIugriz$		13	(47)	0.032	(47)
SN2006aj	Ic-BL	0.0334	-19.33	$UBVri$	u,b,v,uvw1,uvm2,uvw2	22	(26),(8),(24)	0.2	(91)
SN1998bw	Ic-BL	0.0085	-18.93	$UBVRI$		21	(9),(10)	0	(10)
SN2009bb	Ic-BL	0.01	-18.72	$BVugri$		12	(38),(18)	0.540 (0.200)(b)	(88)
SN2007ru	Ic-BL	0.0155	-18.65	$UBVRI$		5	(32),(8),(21)	0	(32)
SN2012ap	Ic-BL	0.0121	-17.66	$BVRI$	u,b,v,uvw1	12	(58),(21)	0.4	(58)
SN2002ap	Ic-BL	0.0022	-17.07	$UBVRI$		23	(13),(14),(8)	0.008	(14)

Equations

$$k(x,x^{\prime })=A\left(1+\frac{\sqrt{3}(x-x^{\prime })}{\sigma }\right)\exp \left(-\frac{\sqrt{3}(x-x^{\prime })}{\sigma }\right),$$

$$f(t)=\{\begin{array}{ll}\alpha (t-t_0{)}^n& \text{if t≥t0}\\ 0& \text{if t<t0},\end{array}$$

$$f(t)=\{\begin{array}{ll}\alpha (t-t_0{)}^n+\mathcal{N}(t_b,{\sigma }_b)(1-e^{t-t_0}),& \text{if t≥t0}.\\ 0,& \text{if t<t0}.\end{array}$$

$${\mu }_{\mathrm{obs}}=m_B+\alpha x_1-\beta \mathcal{C}-M_B,$$

$$f(A_V)\mathrm{d}{A}_V=e^{-0.016A_V}+0.080\mathcal{N}(0,0.801).$$

Full text

††pubyear: 2019††pagerange: Spectrophotometric templates for core collapse supernovae and their application in simulations of time-domain surveys–A.2

1 Introduction

Current and future time-domain optical surveys are expected to discover many thousands of optical transients and supernovae (SNe) on a nightly basis. This number far surpasses the resources available to spectroscopically confirm the nature of each detected event, and thus new approaches must be developed in order to exploit these samples. Techniques to model the composition of transients in future datasets are needed to guide and prioritise the limited follow-up resources, methods to classify transients based on their (possibly incomplete) photometric light curves are required to identify events of interest for different science goals, and accurate models of the expected contamination from core collapse SNe are needed for future cosmological samples of type Ia supernovae (SNe Ia).

Each of these tasks requires the ability to accurately model populations of core collapse SNe. For example, SN photometric classification techniques critically rely on templates of core collapse SNe, whether these are classical ‘template fitting’ techniques (e.g., Kuznetsova & Connolly, 2007; Sako et al., 2008, 2011; Rodney & Tonry, 2009; Gong et al., 2010), or whether these are machine learning techniques (Ishida & de Souza, 2012; Karpenka et al., 2013; Lochner et al., 2016; Möller & de Boissière, 2019) that must be trained on representative samples of data, often based on simulations.

A similar argument can be extended to SN Ia photometric cosmological analyses. The potential of cosmological analyses using photometrically-classified SNe Ia is significant and robust techniques will be required to unlock the power of future samples. In current SN Ia photometric cosmological analyses, the effect of contamination is assessed using simulation techniques (Kunz et al., 2007; Hlozek et al., 2012; Kessler & Scolnic, 2017), which are sensitive to how well the contamination can be modelled, rather than to the classification techniques used. The recent analysis of the photometric SN Ia sample from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS 1; Jones et al., 2017, 2018) has shown that simulations of core collapse SNe using currently available template libraries and luminosity functions significantly underestimate the apparent core collapse SN contamination in the data, and a significant tuning of the underlying properties of the core collapse SN population is required.

Historically, SN spectral templates have been constructed by combining together multiple observations of different SNe of the same type into some representative average spectral series of the class, interpolating between data to produce a spectral series sampled perhaps daily. This is particularly applicable to SNe Ia (Leibundgut et al., 1991; Nugent et al., 2002; Hsiao et al., 2007), where the object-to-object diversity is relatively limited and quantifiable, and the concept of an average spectrum is straight forward. The advantage of this approach is that a spectrum can be interpolated on any epoch, from which photometry can be synthesised. However, this method is not as suitable for the far more heterogeneous core collapse SNe, where the spectral diversity is larger and the observed data typically sparser, thus making it more difficult (and perhaps conceptually meaningless) to combine data from different events into a single template.

This is well demonstrated by the complex classification scheme for core collapse SNe that depends on their spectral (e.g., Filippenko, 1997) and sometimes photometric (e.g., Barbon et al., 1979; Patat et al., 1994; Arcavi et al., 2012) properties; for a recent review see Gal-Yam (2016). Type II SNe (SNe II) have clear signatures of hydrogen in their spectra whereas type I SNe do not; type Ib SNe (SNe Ib) have helium absorption lines whereas type Ic (SNe Ic) do not; and neither SNe Ib nor SNe Ic have the strong silicon or sulphur lines usually present in thermonuclear explosions. SNe Ic are often further divided according to whether they have broad lines (-BL) in their spectra, indicating very high velocities in the expanding ejecta. Type IIb SNe (SNe IIb; Filippenko et al., 1993) are transition objects that begin as SNe II, but lose the hydrogen and develop increasingly strong helium lines as they evolve. Type IIn SNe (SNe IIn; Schlegel, 1990) are SNe II that present narrow hydrogen emission lines, interpreted as the result of interactions between SN ejecta with circumstellar material (Chugai, 1990; Smith et al., 2008). Thus, for core collapse SNe, a library of spectral templates, rather than a single average template, is more appropriate.

The first library of core collapse SN templates was developed for the Supernova Photometric Classification Challenge (SNPhotCC; Kessler et al., 2010a; Kessler et al., 2010b), with 41 templates based on multi-band light curves of spectroscopically-confirmed core collapse SNe from the Carnegie Supernova Project (CSP), the Supernova Legacy Survey (SNLS), and the Sloan Digital Sky Survey-II (SDSS-II). In the absence of high-cadence spectral time series for these events, the same spectral energy distribution (SED) time series was used for each event, taken from the popular templates of Peter Nugent111https://c3.lbl.gov/nugent/nugent_templates.html, and matched to the observed photometry of each individual event. Ultraviolet (UV) and near-infrared wavelengths for these templates are poorly constrained, and only observer-frame $U$/$u$-band photometry is included. This original library has been improved with additional events and is also used in the Photometric LSST Astronomical Time-Series Classification Challenge (PLAsTiCC; The PLAsTiCC team et al., 2018).

Since the SNPhotCC challenge, the number of published core collapse SNe has significantly increased, with the release of large samples of spectroscopic (Modjaz et al., 2014; Hicken et al., 2017; Stritzinger et al., 2018a; Shivvers et al., 2019) and photometric (Arcavi et al., 2012; Taddia et al., 2013; Bianco et al., 2014; Hicken et al., 2017; Gutiérrez et al., 2017a, b; Stritzinger et al., 2018a; Taddia et al., 2018) data, alongside many single object studies of interesting or unusual events (e.g., Pignata et al., 2011; Terreran et al., 2017; Arcavi et al., 2017a; Bersten et al., 2018; Anderson et al., 2018; Izzo et al., 2019). UV coverage of core collapse SNe has also improved, with observations from the Swift satellite using the Ultra-Violet/Optical Telescope (UVOT) instrument (Roming et al., 2005; Bufano et al., 2009; Brown et al., 2014; Brown et al., 2015), and follow-up programs using the Hubble Space Telescope and GALEX (e.g., Gal-Yam et al., 2008; Ben-Ami et al., 2012, 2015). The construction of templates in the UV is essential for simulating transients at higher redshift, where the rest-frame UV is redshifted into the observer-frame optical.

Taking advantage of these new data, in this paper we present a method to construct a SN spectral template based on individual SN events for which well-constrained multi-band light curves and sparse spectroscopic observations have been measured. This technique is data-driven, and in principle can be generalised to any type of transient. Each final spectral template generated by our method consists of a daily sampled spectral time-series, extended into the UV, and optionally corrected for extinction due to dust in the SN host galaxy. Our goal is then to apply this general technique to all suitable core collapse SNe in the literature, providing a template library representing all suitable core collapse SNe. We also provide the general software implemented to build this library, allowing an easy expansion for SNe published in the future.

The paper is laid out as follows. A detailed description of our techniques is presented in Section 2. In Section 3, we then select from the literature a sample of 67 well-observed core collapse SN events for use in constructing our templates, and discuss various challenges in using the data to build the templates. We then illustrate the use of our new library by combining with published estimates of SN rates and luminosity functions (Section 4.1), and simulating core collapse SN contamination in a large photometric SN survey (Section 4.2.1). We present the results and the conclusions drawn from these tests in Section 5. Where relevant, we assume a Hubble constant of $H_{0}=70$ km s*-1* Mpc*-1* and a flat, $\Lambda$CDM universe with a matter density of $\Omega_{\mathrm{M}}=0.3$.

2 Method

In this section, we present the techniques designed to construct time-series spectral templates from photometric and spectroscopic observations of individual SN events. Fig. 1 shows a schematic overview of the technique. (Our method assumes that any photometric data have been corrected for Milky Way and host galaxy extinction at the first order; we discuss this in detail in Appendix A and when applying the code to data in Section 3.)

We first interpolate the observer-frame SN light-curves observed in multiple filters (Section 2.1), and we use this to estimate for each filter the photometric flux at the epochs on which spectral information for the SN is available. The interpolated photometry is then used to ‘flux-calibrate’ the observed spectra by adjusting the overall spectral shape, but retaining the detailed spectral features (Section 2.2). We then use UV photometry from the same SN event to extend the spectral templates at bluer wavelengths (Section 2.2.1). Finally, we smoothly interpolate between the (UV extended) calibrated spectra so that the final time series is sampled daily (Section 2.3). We discuss each step in turn.

2.1 Light curve fits

We explored various techniques to fit or model the SN observer-frame light curves, using different parametric forms. This included testing different models proposed in the literature (Bazin et al., 2009; Kessler et al., 2010b; Karpenka et al., 2013; Sanders et al., 2015; Vacca & Leibundgut, 1996), and exploring different techniques to fit these parametric forms, including the MultiNest package (Feroz & Hobson, 2008; Feroz et al., 2009) and the mpfit package (Markwardt, 2009). However, we found no single functional form that had both a sufficient flexibility to provide good fits across all filters for all core collapse SNe, but yet that could still be reliably constrained by typical SN data. For this reason, we opted instead to use a non-parametric data-driven interpolation technique, Gaussian processes.

2.1.1 Modelling light curves with Gaussian processes

There are now several examples in the literature of the use of Gaussian processes (GPs; Rasmussen, 2006) to interpolate transient light curves or spectral energy distributions (Kim et al., 2013; Inserra et al., 2018; Angus et al., 2018; Saunders et al., 2018). A GP can be seen as the generalisation of a multivariate Gaussian distribution. Formally, a GP generates data located throughout some domain (time, in the case of light curves) such that any finite subset of points in that domain follows a multivariate Gaussian distribution. As a multivariate Gaussian distribution is fully specified by its mean and covariance matrix, a GP is defined by a mean function and the covariance function, $k(x,x^{\prime})$ (also called the kernel), that relates each point (of the light curve in this case) to any other. The hyperparameters needed to define the kernel are usually determined using a likelihood maximisation routine.

Regression using GPs has several advantages. The first is that GP regression is a non-parametric model, and therefore has the flexibility to interpolate the wide heterogeneity of light-curve shapes inherent in core collapse SNe. Secondly, it is straight forward to include in the model the uncertainties of the observed photometric data points and, consequently, to estimate uncertainties on the interpolated light curve. However, while the ‘function-agnostic’ nature of GP is a strength, the downside is that GPs are also physics-agnostic. As a consequence, negative fluxes, diverging and infinitely brightening light curves are, in principle, allowed. However, these non-physical behaviours typically occur when light curves are extrapolated outside the time-range covered by the photometric data, which can be avoided by restricting to interpolation only.

We use the python package george (Ambikasaran et al., 2015) to perform the GP regression on the light curves in flux space. We use a Matern 3/2 kernel function ($k$) of the form:

[TABLE]

where $A$ is an amplitude factor, $\sigma$ regulates the scale at which correlations are significant, and $x$ is the domain over which the regression is being performed, which in this case is time. Generally we optimise the free hyperparameters of the model (amplitude and scale of the kernel function) by minimising the log-likelihood of the model given the data, and we set the prior on the mean function to a constant zero function.

2.1.2 Early and late phases

The GP interpolation works well at most light-curve phases, but requires extra attention at the start and end of the light curve (where photometric sampling is sparser or absent) to ensure that the interpolation is robust.

At late phases, the SN luminosity is powered by 56Co decay. Sparse photometric measurements are sufficient to constrain the SN evolution in this phase, and so SNe are typically observed less frequently. However, when the data become sparser, the GP interpolation becomes less informative, and with larger uncertainties. We therefore visually inspect each light curve and determine at which phase it becomes dominated by 56Co decay. We then interpolate/extrapolate additional photometric points and propagate their uncertainties from the linear fit for use in the GP interpolation. This ‘oversampling’ of the light curve significantly improves the GP interpolation, without adjusting the original data around peak.

At early phases we use a similar approach and over-sample the rising part of the light curve, filling the gap between the estimated explosion day and the first photometric data point. We fit the early data flux $f$ as a function of time $t$ with a widely used parametrization of

[TABLE]

where $\alpha$ is a normalisation coefficient, $t_{0}$ is the time of explosion, and $n$ is the power index of the rising light curve.

For most of the SNe, an accurate estimate of $t_{0}$ can be found in the literature, typically derived from light-curve modelling or from the analysis of early spectroscopic observations (e.g., Drout et al., 2016). Thus, we generally fix the parameter $t_{0}$. When the explosion epoch is uncertain, we treat $t_{0}$ as a free global parameter and use non-detections and the SN date of discovery as lower and upper bounds on $t_{0}$. The power index $n$ is fixed only when the available early photometry is too poor to constrain the fit. For stripped envelope SNe, we assume a power law of the form $f(t)\propto t^{1.5}$, expected for the cooling of the shock-heated expanding ejecta (Piro & Nakar, 2013). For hydrogen-rich SNe, we assume$f(t)\propto t^{0.935}$, (González-Gaitán et al., 2015).

For some SNe considered (SN 1987A, SN 1993J, SN 2006aj, SN 2013df, SN 2011dh, SN 2011fu) an initial shock-breakout is detected in the light curves, presenting as an initial peak in the photometry. In such cases, the fitted functional form used is:

[TABLE]

where $t_{b}$ and $\sigma_{b}$ are the peak and width of the initial bump in the light curve.

Fig. 2 shows an example of a multi-colour light curve interpolation using GPs for the SN iPTF13bvn (Fremling et al., 2014). Equivalent figures for all other events in our sample are available as online-only figures.

2.2 Generating spectrophotometric data

Observed SN spectra are rarely spectrophotometric, no matter how carefully the observations and reductions are performed. Differential slit losses during the observation, centring errors of the object on the slit, or non-photometric observing conditions not only affect the overall normalisation of the spectrum, but can also introduce a smooth, wavelength-dependent distortion of the continuum (Buton et al., 2013).

Multi-band photometry allows a correction to be made for many of these effects: it can set the overall flux normalisation, and also provides information on the broad-band colours of the spectrum. We therefore adjust the observed spectra so that they match, in colour space, the photometry interpolated at the epoch of spectral observation. This correction is sometimes referred to as ‘mangling’ the spectrum (e.g., Hsiao et al., 2007; Conley et al., 2008), and consists of determining a smooth wavelength-dependent function that, when multiplied by the original spectrum, produces a spectrum with the correct colours and normalisation. In order to calculate this smooth wavelength-dependent function we again use GPs.

For each SN spectrum, we first calculate the observed flux in each of the filters for which a photometric light curve is available. We compare these synthetic fluxes with the photometry interpolated from the light curve (Section 2.1), and determine the ratio of the two, propagating the uncertainties from the GP light curve interpolation and from the spectrum. We place each calculated ratio at the filter effective wavelength of the original spectrum, and use GPs to interpolate as a function of wavelength, obtaining a smooth wavelength-dependent calibration function. We use again a Matern 3/2 kernel (Eq. 1), fixing the scale to 300 Å.

We do not use the common approach of using spline interpolation in the mangling. We explored this approach in detail, but found it is difficult to account correctly for the observational uncertainties in the photometry. The spline interpolation can also produce discontinuities in the mangling function for SNe for which photometry in neighbouring filters with close effective wavelengths is available (e.g., SNe which have photometry in both SDSS $r$ and Bessell $R$). Fig. 3 shows an example of a mangled spectrum and the interpolated mangling function.

2.2.1 Extending the wavelength coverage of the spectra

During the mangling procedure described above, any near-UV or near-infrared (near-IR) data from the SN are not used, as the spectroscopic data typically cover the optical wavelength range222Less then 10 per cent of the spectra considered in Section 3 have coverage below 4000 Å and only 50 per cent cover $>$9000 Å.. However, near-UV coverage in the SN rest frame is critical to simulate SN events at higher redshift (i.e., typically at $z>0.2$). We therefore use near-UV photometry observed for the same SN to smoothly extrapolate the spectra to lower wavelengths.

We experimented with smoothly extrapolating the spectrophotometric data by fitting a black body function to the optical and near-UV photometry. While this works for the rather featureless very early spectra of some SN subtypes, it provides a poor near-UV representation for most of the SN events we consider. The complex absorption due to ionised metals present in the SN envelope typically dominates the UV part of the SN SED during most of its evolution, and its modelling is non trivial (Gal-Yam et al., 2008; Bufano et al., 2009; Ben-Ami et al., 2012, 2015). We therefore decided to use a less rigid approach and again implement GPs. In this case, we combine photometry (i.e., integrated flux versus time, $t$) and spectrophotometry (i.e., flux density versus wavelength, $\lambda$) into a two-dimensional grid on which we can interpolate a flux surface $f(t,\lambda)$.

The implementation of a two-dimensional interpolation of this flux surface $f(t,\lambda)$ presents several advantages. First, it guarantees that the near-UV extension is continuous in time; extending each spectrum individually can result in large discontinuities in the final spectral time series in the near-UV. In addition, it is common to have spectroscopic observations of the same SN from different facilities. As a result, the wavelength coverage of the available spectra can change significantly from spectrum to spectrum, with some spectra extending to $\sim 3200$ Å, while other ending at 4000-4500Å. A two-dimensional interpolation allows us to perform an extension including any information from adjacent spectra, and provides a final template that is continuous in time. Moreover, a two-dimensional model not only allows us to extend the spectrophotometry into the near-UV or near-IR, but also to simultaneously smoothly interpolate between the spectra and extrapolate additional spectra when the available data are too sparse (see Section 2.3).

In order to interpolate the flux surface $f(t,\lambda)$, both photometric and spectrophotometric data are combined and smoothed on a grid with a 60 Å binning in wavelength (to uniformly sample all the calibrated spectra and reduce the number of data points used to train the GP model). An example is shown in Fig. 4. The UV integrated flux is included and placed on the wavelength axis at the mean wavelength of the filter. The kernel used for GP interpolation is a two-dimensional Matern 3/2 Kernel function (Eq. 1), where $x$ is now the two-dimensional vector $\mathbf{x}=\{x_{t},x_{\lambda}\}$. The GP hyper-parameters, $\sigma_{\lambda}$ and $\sigma_{t}$ are fixed to 100Å and 30 days respectively (optimising these hyper-parameters does not affect the results significantly and is computationally expensive).

The two-dimensional GP interpolation requires a prior for the mean function to prevent non-physical behaviour (i.e., negative fluxes or a GP model rapidly decreasing to zero when the data coverage is poor) from occurring. Building this prior is not trivial given the diversity of core collapse SNe in terms of brightness, shape and colours. In order to do this we use observed data from the 67 core collapse SNe that we will introduce in Section 3, examining their colour evolution in the optical and near-UV after dust extinction corrections (these SNe are all at $z<0.03$ where $k$-corrections are small). To first order, SNe II and fast rising SNe IIn (hydrogen-rich SNe) present a similar colour evolution, despite having different absolute brightnesses. The same is true for the stripped envelope SN classes of SNe IIb, SNe Ib and SNe Ic/Ic-BL (see also Drout et al., 2011; Stritzinger et al., 2018b) and for the sub-class of slow rising SNe IIn like SN 2006aa and SN 2011ht.

First, we measure, for each SN and for each optical and near-UV filter $X$, the evolution of $F_{X}(t)/F_{V}(t)$, where $F_{X}$ and $F_{V}$ are the GP interpolated light-curves in filters $X$ and the $V$-band respectively. We then average $F_{X}(t)/F_{V}(t)$ over all SNe within the same sub-class (hydrogen rich, stripped envelope or slow SNe IIn). We also measure, in each filter, the average wavelength weighted by a typical hydrogen-rich, stripped-envelope or slow Type IIn SED. We interpret each average colour evolution $\langle F_{X}(t)/F_{V}(t)\rangle$ as a monochromatic measurement of colour at the respective average wavelength, and we smoothly reconstruct the ‘colour surface’ $\langle F_{X-V}(t,\lambda)\rangle$. For each SN, the prior used to interpolate the flux surface $f(t,\lambda)$ is calculated multiplying the ‘colour surface’ $\langle F_{X-V}(t,\lambda)\rangle$ by the $V$-band light curve of the SN itself. This means that for each SN we build a different prior that has the colour properties of the sub class of SNe it belongs to, but is still normalised to the apparent brightness of the SN considered.

We test the effects of our prior using some of the best observed SNe (e.g., SN 2013by and iPTF13bvn). We perform the two-dimensional GP interpolation with and without the prior, using the full set of data available for the SN, and then removing some information (i.e., simulating reduced UV coverage or sparser spectroscopy). We find that the choice of the prior improves the results when sparser data or less UV coverage is available, without significantly overestimating or underestimating the final near-UV extension, and without affecting the accuracy of the spectrophotometry.

Once the prior is built, the flux surface $f(t,\lambda)$ is interpolated and it is then straight forward to extrapolate for each optical spectrum the near-UV extension. Fig. 5 shows an example of spectra observed for the SN Ib iPTF13bvn and extended applying the technique presented in this Section. In Fig. 4 we show, for the same SN, a graphical representation of the flux surface $f(t,\lambda)$ as reconstructed from the two-dimensional GP.

2.2.2 Refining the calibration of the extended spectra

The extended spectra reconstructed from the two-dimensional interpolation provide a good approximation of the underlying near-UV SED. However, when building the two-dimensional surface, we interpret the photometric measurements as monochromatic flux density measurements at the filter mean wavelength. Ideally, the effective wavelength – $\lambda_{\mathrm{eff}}$, the average wavelength of the filter weighted by the SED – should be used instead, but this is not known a priori as the SED is not known. This can lead to inaccuracies, particularly in the near-UV when using photometry from Swift (Brown et al., 2016), where the filter transmission functions present significant tails into the optical. This can be relevant for some core collapse SNe with strong UV variability or at very early phases. As a consequence, the (time dependent) filter effective wavelength can be significantly different from the filter mean wavelength as the SED changes, and as a result, fluxes synthesised from our extended spectra in the Swift filters can be up to 40-50 per cent discrepant when compared to the observed photometric fluxes.

We address this by adjusting the near-UV fluxes to remove the optical contribution from the filter red leaks before including them in the two-dimensional grid for interpolation. We do this by integrating the spectra over each filter’s optical tail, subtracting from the original Swift photometry, and then using this subtracted photometry in the two-dimensional grid for the interpolation. After the spectral time series is extended, we then apply a further correction to the final extended spectra. For each extended spectrum, we compute synthetic photometry in the near-UV and optical filters, and compare with the original photometry observed in the same filters. We then adjust the extended spectra until the synthetic and observed photometry match. We apply this process iteratively as $\lambda_{\mathrm{eff}}$ for each filter changes as the spectrum is mangled. We refer to this additional correction as ‘remangling’. At the end of this process, the extended spectra are robustly flux-calibrated at all wavelengths.

During the remangling, the GP model provides not only a prediction of the mangling function, but also the relative full covariance matrix. This includes the uncertainties on the input photometry (added in quadrature to the diagonal of the covariance matrix of the GP model) and provides a robust estimate of the (highly correlated) uncertainties on the final flux calibrated spectrum.

2.3 Final templates

The most widely used SN simulation packages usually simulate SN light curves by synthesising broadband photometry from a spectral time series model, and then linearly interpolating between this synthetic photometry to generate a light curve at any required cadence. For this reason, the optimal time sampling for our final time series SED templates is approximately one or two days, for a total phase coverage from a few days prior to maximum brightness to around a hundred days after maximum brightness. The observed photometry for our SN events typically meet this quality of sampling. However, the spectroscopic observations generally have a significantly lower cadence and poorer phase coverage.

We therefore increase the time sampling of our templates by smoothly interpolating between the calibrated spectra. This smooth interpolation is performed as described in Section 2.2.1, and an example is shown in Fig. 4. Additional spectra are extrapolated from the interpolated flux surface and then remangled as described in Section 2.2.1, ensuring that the colours and normalisation of the synthetic extrapolated spectra data are accurate.

As a final step, we transform each template to the rest-frame using the measured heliocentric redshift. When necessary, corrections for peculiar velocities are applied. Each SN template is daily sampled and its phase is calibrated arbitrarily by defining day zero as the day at which the synthesised $BVRI$ pseudo-bolometric light curve reaches its maximum value.

The result is a daily sampled and de-reddened, rest-frame spectrophotometric template. These templates can be used for simulations of SN light curves up to $z\approx 1$, in any optical filter system and assuming any Milky Way or host galaxy dust reddening. We also provide templates for which host-galaxy extinction has not been corrected. In Fig. 6, we present synthetic photometry simulated at different redshifts in the LSST filter system using the SN template built from the SN Ib iPTF13bvn.

These template SEDs can be easily integrated in SN simulation packages such as SuperNova ANAlysis software (snana; Kessler et al., 2009) or sncosmo (Barbary et al., 2016). The Python package used to generate the templates is also available333https://github.com/maria-vincenzi/PyCoCo, along with examples and tutorials.

3 Data

We now describe the data from which we can use the techniques described in the previous section to construct a series of core collapse spectral templates. We do this using 67 core collapse SNe for which high-quality photometry and spectroscopy has been published. We describe the selection criteria we applied to build this dataset, and various other data preparation procedures.

3.1 Description of the data set

3.1.1 Selection criteria

We select core collapse SNe from the literature according to the following criteria:

• •

The published SNe must have photometric data in at least the three optical filters of $B$,$V$ and $R$ (or alternatively $g$,$r$ and $i$), with a well-constrained light curve interpolation and at least one epoch prior to peak brightness in the $V$ (or $g$) band.

• •

The SNe must have photometry in at least one near-UV filter, i.e., a filter that samples the wavelength range below 4000Å in the SN rest-frame. This near-UV photometry must have at least one epoch within two days of the $V$/$g$ band peak.

• •

The SNe must have at least five spectra, with at least one observation within three days of the $V$/$g$-band peak, and at least one observation after +15 days.

We do not apply any redshift constraints, and in principle our technique can be applied to SNe observed at any redshift. However, multiple spectroscopic follow-up observations of SNe become expensive at redshifts above $z\sim 0.05$, and thus mostly nearby SN events have been included in this work.

We found 67 core collapse SN events in the literature that satisfy these quality cuts. This sample is presented in Table 2. Where available, we download the spectra from the WISeREP repository444https://wiserep.weizmann.ac.il/ (Yaron & Gal-Yam, 2012) and the photometry from the ‘Open Supernova Catalog’ (Guillochon et al., 2017).

3.1.2 UV data

Since all SNe selected are at low redshift ($z<0.03$; Table 2), near-UV photometry in the SN rest-frame corresponds to observations in the observer-frame $U$ band or bluer filters. For 36 of the 67 SNe selected in this work, UV photometry is measured by Swift/UVOT (Brown et al., 2014) using four UV filters (and two optical filters) for imaging, and two grisms for low-resolution spectroscopy. The central wavelengths for the near-UV filters are 3465Å, 2600Å, 2246Å and 1928Å for the filters $u$, $uvw1$, $uvm2$, $uvw2$ respectively.

3.1.3 Host galaxy extinction estimates

Providing a library of host galaxy extinction-corrected core collapse SN templates is a significant advantage, especially when using these templates in simulations, as arbitrary amounts of extinction can then be added. Thus the input to our template construction method in Section 2 should be extinction-corrected photometry.

However, estimating SN extinction due to local host galaxy dust is challenging as indicators of dust extinction are difficult to measure. We discuss various methods and resources used to estimate the SN host galaxy extinction in Appendix A. In Table 2 we report for each SN the values adopted for reddening due to dust in the SN host galaxy, $E(B-V)^{\mathrm{host}}$. We note that our code will make the required corrections to the SN photometry provided the extinction estimates (i.e., $E(B-V)$) are provided.

3.1.4 Supernova classification scheme

In the application of our library of templates for simulations, and for training supervised machine learning classifiers, we require a broad classification scheme. Given the limited number of objects available in our dataset, we group our templates into six spectral sub-types: SN Ib (13 SNe), SN Ic (7 SNe), SN Ic-BL (6 SNe), SN IIb (11 SNe), SN II (23 SNe), SN IIn (7 SNe) and SN 1987A.

We do not apply the historical distinction between SNe IIP and SNe IIL, as recent analyses show reduced evidence for such a separation (Anderson et al., 2014). The sample of 23 SNe II considered in this work present a continuum of decline rates, from small decline rates and almost flat light curves, to fast-declining events fading at around one magnitude per 30 days.

We also include in our sample more recently identified sub-classes of peculiar SNe, such as SNe Ibn (SN 2010al) and other peculiar SN IIn events (IIn-pec) such as SN 2009ip. SNe Ibn/IIn-pec represent very interesting classes of rare transients: their luminosities can be comparable to SNe Ia and their light curves are not necessarily as wide as normal SNe II. However, the lack of unbiased samples of these classes of transients, and their diversity, makes it difficult to measure their global properties (e.g., their luminosity distribution or relative rates), and therefore, we combine SNe Ibn and SNe IIn-pec into a single group of SNe IIn. In Table 2, we report the list of SNe included in this work and the spectroscopic type.

3.2 Implementation

The data that we use to build the spectral templates are generally used as published. However, in the implementation phase it is sometimes necessary to remove observational and astrophysical signatures that do not have a SN origin, or that are not relevant for the purpose of these templates. We also apply some basic data conversions and noise reduction procedures as follows:

As the light curve interpolation is always performed in flux space, we convert photometric measurements from magnitudes to flux densities. For the conversion of photometry to flux densities we use our own custom-written software555https://github.com/chrisfrohmaier/what_the_flux. 2. 2.

We correct the photometry for Milky Way and host galaxy extinction using a Cardelli et al. (1989) dust law with $R_{V}=3.1$. For Milky Way extinction we use the Schlegel et al. (1998) dust maps. We discuss in detail our choice of $R_{V}^{host}=3.1$ in Appendix A. 3. 3.

Some older events (e.g., SN 1987A, SN 1993J) do not have published uncertainties on their magnitude measurements. We arbitrarily estimate the uncertainties as 10 per cent of the fluxes. 4. 4.

Photometric observations taken on the same night and in the same filter are averaged in order to reduce uncertainties in the individual measurements. 5. 5.

Using a sigma-clipping algorithm, we remove spectral features clearly associated with a host galaxy (e.g., nebular emission lines) or due to the atmosphere (telluric features). 6. 6.

We smooth each spectrum by applying a Savitzky-Golay filter (Savitzky & Golay, 1964) using a window of 100 Å.

Once the data has been pre-processed as described, we apply on them the technique presented in Section 2. In around 6 per cent of the light curves considered in our data set, the signal-to-noise of the photometry is sufficiently low such that the GPs over-smooth the data. In such cases, we manually fix the scale of the kernel or use a non-zero prior on the mean function when performing the light-curve fitting described in Section 2.1.

After testing our technique on the sample of 67 SNe described in Section 3.1, we compare the published observed photometry with the photometry synthesised from the final flux-calibrated, daily-sampled and near-UV extended spectral time-series. Our results are presented in Fig. 7. On average the observed photometry is recovered to within 0.02 magnitudes in most of the optical filters and to within 0.1 magnitudes in the near-UV.

4 Application to simulations of SN surveys

Our new library of core collapse SN templates has many potential uses in photometric SN classification and the simulations of optical transient surveys. In this section, we demonstrate their application to SN simulations by integrating our library of core collapse SN templates into the SN light-curve simulation package snana. This requires two further key ingredients in simulating SNe, no matter what the scientific goal: the luminosity function used to set the overall brightness distribution of the simulated events, and the relative rates that different SN subtypes are simulated with. We then present an example application: the simulation of contamination by core collapse SNe in cosmological surveys of SNe Ia.

4.1 Luminosity functions and relative rates

The luminosity function (LF) describes the intrinsic distribution of brightnesses (absolute magnitudes) of SNe. Due to selection biases inherent in all astronomical surveys (or compilations/catalogues), it is difficult to measure accurately, with two main core collapse SN luminosity functions published in the last decade: Li et al. (2011) and Richardson et al. (2014).

The brightness distributions presented in Li et al. (2011, hereafter L11) are measured from a volume-limited sample of 92 core collapse SNe from the Lick Observatory Supernova Search (LOSS; Leaman et al., 2011). They are constructed in the $R$-band filter, and calibrated to the Landolt (1992) system from mostly unfiltered images with a precision of about five per cent (Li et al., 2003). No $k$-corrections were made as the SNe all lie within 60 Mpc and thus $k$-corrections are small. In order to correct for potential biases in the volume-limited sample, L11 estimated for each SN the completeness of the LOSS survey within the fixed cutoff distance, given the brightness of the event. The peak absolute magnitude and light-curve shape of each SN are recovered by performing template fitting of the $R$-band light curve then used in the computation of the completeness (see also Leaman et al., 2011). Different templates are used for SNe of different types. Each SN in the LF sample is then weighted by the reciprocal of this completeness.

Recently, Shivvers et al. (2017) revised the classifications of the SNe in the LOSS sample, and presented an updated measurement of the relative rate of each SN sub-type. In particular, SNe IIL and IIP are now grouped into a single population (mirroring our templates), and calcium-rich SNe (Perets et al., 2010) and SNe Ic-BL are considered distinct classes. The classification of 15 SNe has also been corrected using additional spectroscopic observations. We use this re-classified LOSS sample to remeasure the LF of each core collapse SN sub-type, with SNe Ic, Ic-pec and Ic-BL grouped into the same sub-class due to the small numbers of events included in these classes (the relative fraction can be preserved in any simulations). We also exclude two SN IIn now reclassified as SN imposters.

In principle, for each re-classified SN a new template fit and the re-computation of the completeness following L11 would be required. However, L11 used the same template for SNe Ib and Ic, and thus SNe Ic re-classified as SNe Ib do not require new completeness calculations.

Other SNe that do require a new light curve fit and for which the completeness should be recalculated are partially or completely excluded. This is the case of 5 SNe. SN 2003bw, SN 2003br and SN 2005mg are re-classified from SNe II to SNe II/IIb. We can include them in SNe II LF but not in SNe IIb LF. SN 2004C and SN 2005lr are re-classified from Ic to IIb, therefore we can not include them in the SN IIb LF. These approximations are likely to bias the LF for SNe IIb. However, an accurate estimate of LFs goes beyond the purpose of this work. A comparison between the original L11 LFs and revised ones is shown in Table 1.

We also consider the LFs of Richardson et al. (2014, hereafter R14). R14 measured LFs from a larger sample of 211 core collapse SNe, combining the Asiago SN Catalog (Richardson et al., 2002, 2014) and supplementary magnitude-limited SN samples. These LFs are presented in the $B$ band, and have additionally been corrected for host-galaxy extinction. However, only half of the SNe considered (113 of 211) have specific information about host galaxy extinction, with an average $B$-band extinction of 0.292 mag assumed for the other events. Intrinsic biases in the SN sample are corrected using a bootstrap technique, and LFs measured from both the bias-corrected SN sample and from a volume-limited but not bias-corrected SN sample. We test in our simulations the bias-corrected LFs (also reported in Table 1).

To determine the relative fraction of each core collapse sub-type, we use the revised nearby SN relative rates presented in Shivvers et al. (2017, see also Table 1). These relative rates are assumed across all our simulations.

The LFs and the relative rates used in this section are all based on local SN measurements. In our simulations, we assume these LFs and relative rates also apply at higher redshift, i.e., we neglect any redshift evolution in these quantities. The only redshift dependency we assume is that of the global volumetric rate of all core collapse SNe.

4.2 Simulating supernovae in cosmological surveys

SNe Ia can be used as standardisable candles to estimate distances in the universe (Riess et al., 1998; Perlmutter et al., 1999). Historically, samples of SNe Ia have been constructed where every SN Ia event has been both spectroscopically confirmed and has a spectroscopic redshift measurement (e.g., Betoule et al., 2014; Scolnic et al., 2018). However, contemporary and future surveys will be unable to spectroscopically confirm every SN Ia event detected due to a lack of follow-up resources, and will therefore be reliant on photometric classification techniques based on a multicolour light curve and possibly a spectroscopic redshift of the SN host galaxy.

Contamination from events that are not SNe Ia is therefore potentially a serious problem (Campbell et al., 2013; Jones et al., 2017). Simulations can be used to estimate selection effects and the level of this contamination (Kessler & Scolnic, 2017); a recent example can be found in Jones et al. (2017, hereafter J17) and Jones et al. (2018) using the the Pan-STARRS 1 (PS1) photometric SN sample.

These simulations have analysed the ‘Hubble residual’ distribution for the SNe Ia in the sample. The Hubble residual is defined as $\mu_{\mathrm{obs}}-\mu_{\mathrm{model}}$, where $\mu_{\mathrm{obs}}$ is the observed distance modulus for each SN event, and $\mu_{\mathrm{model}}$ is the theoretical distance modulus expected for the redshift of each event given a set of cosmological parameters. $\mu_{\mathrm{obs}}$ is defined as

[TABLE]

where $m_{B}$ is the rest-frame $B$-band apparent magnitude at maximum light (or more accurately the log of the light-curve amplitude in the light-curve fit), $x_{1}$ is the SN ‘stretch’ parameter, $\mathcal{C}$ is the SN colour estimator, and $\alpha$, $\beta$ and $M_{B}$ are parameters defining the stretch–luminosity and colour–luminosity relationships for the SN Ia population being studied (e.g., Tripp, 1998; Astier et al., 2006). The $x_{1}$ and $\mathcal{C}$ estimators are typically determined using the SALT2 light-curve fitter (Guy et al., 2007).

The set of simulations implemented for the cosmological analysis of the PS1 sample is presented in detail in J17. The parent distributions of $x_{1}$ and $\mathcal{C}$ used to simulate SNe Ia are determined by applying the method of Scolnic & Kessler (2016) to the PS1 spectroscopic SN sample (Rest et al., 2014) with additional adjustments, taking into account additional selection biases on the photometric sample. Core collapse SNe are simulated using the set of 41 templates from the SNPhotCC library, with the addition of four new templates of SNe IIb built using more recent observations (models for this SN sub-type were not included in the original SNPhotCC library). Templates of sub-luminous SN 1991bg-like SNe Ia are also included. This library of templates is not corrected for host extinction and is matched to the original LFs taken from Li et al. (2011). This approach of modelling the core collapse SN population assumes that the reddening distribution in the data, the reddening distribution in the templates, and the reddening distribution in the sample of SNe used to measure the LFs are equal and representative.

A significant discrepancy is observed between the Hubble residual distribution measured from the SN data, and that predicted from the simulations: a third of the SNe observed in the data with Hubble residuals in the range $0.5<\mu_{\mathrm{obs}}-\mu_{\mathrm{model}}<1.5$ mag are not present in the simulations. In J17, this is accounted for by adjusting the core collapse LFs, brightening the Li et al. (2011) values by $\simeq 1$ mag, and fine-tuning the width of the LF (these adjusted LFs are also reported in Table 1). As discussed in J17, these adjustments are made after analysing the peak absolute brightnesses of the PS1 light curves identified as SNe Ib/Ic or SNe II by the photometric classifier psnid (Sako et al., 2011), after specific SALT2 cuts are applied to the sample, and are not necessarily representative of intrinsic SN LFs.

4.2.1 Implementation of the simulations

We repeat the J17 simulations of the PS1 photometric SN sample using our new template library to demonstrate its applicability to simulations of SN surveys. Our simulations use the same snana configuration files presented in J17, with the exception of revising the approach used to model and reproduce the core collapse population. We test three alternative ways to generate core collapse SNe:

Using our uncorrected templates, matched with the revised L11 LFs (hereafter referred to as V19+L11), 2. 2.

Using our host-extinction corrected templates matched with the R14 LFs, simulating additional host extinction (‘V19+R14 LFs’), 3. 3.

Using our uncorrected templates, matched with the LFs adjusted by J17 (‘V19+J17’).

For each approach we generate 25 independent simulations, and our results for each case are then an average of these 25.

The method used to match the templates to the relative LF is straight forward. For each of the SN sub-types considered in our classification scheme, we compare the (biased) luminosity distribution measured from our set of templates to the corresponding LF. Sub-type dependent magnitude offsets and smearing factors are applied to the templates so that the peak brightnesses of simulated SN events have the correct distribution. By applying sub-type dependent magnitude offsets, we correct the overall brightness of each sub-class of templates, thus retaining the relative brightness between individual templates of the same sub-type. As a result, correlations intrinsic to sub-classes of templates (e.g., faster declining SNe II are on average brighter than slower declining SNe II) are preserved in the simulations.

The relative rate of each template is calculated so that the core collapse SN fractions presented in Shivvers et al. (2017) are correctly reproduced in the simulations. Within a sub-type, each individual template is considered to be equally probable. For simulations using templates corrected for host galaxy extinction, the effect of host galaxy dust is simulated using the dust extinction distribution presented in Hatano, Branch & Deaton (1998, see Appendix A).

We compare the results from these three approaches with the core collapse contamination modelling adopted in J17 after the fine-tuning of the LFs (hereafter referred as J17 templates with J17 adjusted LFs). This gives four families of simulations in total, three of which we run ourselves. The input files used to generate each simulation can be found in the package snana.

4.3 Results and Discussion

We now turn to the results from our simulations. For each of the simulations, detected SNe are fitted using the SALT2.4 model (Guy et al., 2010), and the SALT2-based cuts of Betoule et al. (2014) and J17 are then applied. In J17, over 25 independent simulations, the average number of simulated SNe Ia passing SALT2-based cuts is $1021\pm 30$, which correspond to $88.4\pm 3.7$ per cent of the number of SNe observed in the PS1 sample after SALT2-based cuts (1153 SNe $\pm$ a Poisson error of 34). Assuming that simulations of SNe Ia are correct and that each event in the data is either a SN Ia or a core collapse SN, we expect to reproduce with simulations a fraction of core collapse of $11.6\pm 3.7$ per cent.

In Fig. 8 we present for each simulation the predicted core collapse contamination and the distribution of Hubble residuals. (In the discussion that follows, we remind the reader that the Hubble residual is defined as $\mu_{\mathrm{obs}}-\mu_{\mathrm{model}}$; i.e. brighter SNe have a negative Hubble residual.)

Comparing our results to the PS1 data we observe that:

• •

Simulations performed using revised L11 LFs and templates not corrected for host extinction (V19+L11) predict a core collapse contamination of $3.7\pm 0.5$ per cent. The number of SNe in the fainter tail of the Hubble residuals is underestimated by a factor of approximately two. This reproduces the result of J17 before the fine-tuning of the original L11 LFs, and shows that L11 LFs, used either in their original or revised form, underestimate the number of observed contaminants.

• •

Simulations performed using R14 LFs and de-reddened templates not only underestimate the number of SNe in the fainter tail of the Hubble residuals but also overestimate the number of bright events (Hubble residuals with $<-1$ mag) by a factor of three. The overall predicted contamination is $9.5\pm 0.9$ per cent. The high fraction of bright events suggest that either the R14 LFs are over-represented in bright objects, or the dust extinction applied in the simulations is underestimated (globally or for particular SN sub-types).

• •

Simulations using the LFs adjusted by J17 (V19+J17 and J17 simulations) better reproduce the overall Hubble residual distribution as expected, and predict a contamination of $7.5\pm 0.8$ and $7.7\pm 0.9$ per cent respectively. However, we note that when our library of templates is used instead of the original J17 library, the agreement between the simulations and data slightly decrease. This suggests that LFs fine-tuned by J17 should be used with caution as they are fine-tuned to a specific set of core collapse templates.

In Fig. 9 we compare simulated and observed Hubble residual distributions in different redshift bins ($z<0.3$, $0.3<z<0.4$ and $z>0.4$). We present separately the effects of the different LFs while fixing the template library (upper plot in Fig. 9) and then compare simulations implemented using the same LFs (J17 adjusted) but different templates libraries (lower plot in Fig. 9). Different redshift ranges present different issues:

• –

At lower redshifts ($z<0.3$) none of the simulations generated using our templates agree with the data. Despite the fact that the LFs tested are significantly different in terms of brightness (see Table 1), the number of SNe for Hubble residuals $>0.5$ mag is systematically under-predicted by a factor of two to three (corresponding to 20-25 ‘missing’ SNe). Even using de-reddened templates and simulating a range of dust extinctions (V19+R14 set of simulations), we do not reproduce the observed contamination despite increasing significantly the colour diversity of the simulated light-curves. At low redshifts, fainter SN events (e.g., stripped-envelope SNe or extincted SNe) are more likely to be detected; therefore the diversity of the template library, the fraction of peculiar events included, and the extinction distribution assumed for each SN sub-type, all have a significant impact on the final results.

• –

At higher redshift the fraction of detected core collapse SNe is reduced and the agreement between the data and simulations improves significantly. At $z>0.4$, we note that the expected core collapse fraction is larger than $12$ per cent (with $366\pm 16$ simulated SNe Ia in this redshift bin, and 417 observed SNe in the PS1 sample). This raises some concerns about the modelling of SNe Ia in this redshift bin. We note that simulations using the R14 LFs can predict such a high contamination, but significantly over-predict the number of bright SNe (Hubble residuals of $<-1$ mag). Simulations using the J17 LFs better reproduce the Hubble residuals in this redshift bin, but here we note that when our ‘near-UV extended’ templates are used, the predicted core collapse contamination is two times larger than that produced with other templates. This suggests that the near-UV extension is important for simulations of core collapse contamination at higher redshift.

Finally, Fig. 10 shows the contribution of each core collapse SN type to the final simulated Hubble residual distributions. Several interesting trends emerge, highlighting the importance of the template library used in the simulations. For example, the number of SNe Ib and SNe Ic with Hubble residuals between $-1$ and $+1$ mag is larger for simulations in which the R14 LFs and de-reddened templates are used. Two effects drives this result. Firstly, a significant fraction of the SNe Ib and Ic entering our templates are highly reddened ($E(B-V)>0.2$; see Fig. 14), and many simulated SNe Ib and Ic do not pass the SALT2 colour cuts if templates that include host extinction are used. Second, the R14 LFs for SNe Ib ans SNe Ic are brighter than that of L11 or J17, even accounting for the differences in filters in which the two functions are constructed ($B$ band versus $R$ band).

We also note that more SNe IIn are detected in simulations using our templates, even when the same LF is used (third and fourth panel in Fig. 10). This is due to the larger diversity of SNe IIn captured in our library of templates, with both fast- and slow-declining SN IIn events included. In simulations where the L11 revised LFs and R14 LFs are used, a larger number of SNe IIn is expected as the intrinsic brightness of SNe IIn is larger than in the adjusted LFs of J17 (see Table 1). Finally, we observe that the Hubble residual distribution of SNe II from simulations using our library of templates is shifted by approximately 0.5 magnitudes when compared to simulations using the J17 templates. This is true even when the same underlying LFs are used (third and fourth panel in Fig. 10) and is related to differences in the colour and stretch properties of the templates.

To conclude, analysing the Hubble residual distribution (globally, for different redshift bins and for different SN sub-types) is an important test for cosmological analyses. However, Hubble residuals are (by definition) conflated with SN parameters such as stretch and colour (equation 4), and various cuts based on these SN properties are also applied. Thus, the analysis of Hubble residuals alone cannot provide a fundamental understanding of whether it is the luminosity functions, or the colour and stretch distributions and associated cuts, or both, that are driving the differences between observed and simulated samples. We leave for a future article a detailed comparison between various global properties of observed and simulated SNe, using the samples both before and after SN parameter based cuts. However, understanding what drives the discrepancies between the data and simulations is crucial and provides important guidance in the design of realistic simulations of SNe of multiple types, essential to reduce systematic uncertainties due to core collapse contamination in future SN Ia analyses.

5 Summary and future applications

We have presented a library of 67 spectral time-series templates of core collapse SNe. Compared to existing core-collapse spectral models, our new set of templates represent a significant improvement in several ways. We highlight the following key features of our library:

Spectral templates included in the library are generated from multi-band photometry and only sparsely-sampled spectral data. 2. 2.

Each template is generated from a single SN event, retaining the event-to-event diversity in the final library. No SED template or model is assumed, and the final daily-sampled spectral time series for each SN is reconstructed in a fully data-driven fashion, with a heavy emphasis on Gaussian processes. 3. 3.

Each template is extended in the near-UV to $\simeq$1600Å using near-UV photometry. This is critical to simulate SN light-curves at redshifts greater than 0.2. 4. 4.

Each template has (optionally) been corrected for dust extinction in the SN host galaxy. As a result, different levels of host galaxy extinction, or different dust models, can be applied to the templates in simulations, and a realistic range of reddened SN light-curves can be reproduced.

The techniques and code used to construct the spectral library are open-source, and can, in principle, be applied to any type of transient with well-constrained multi-band photometry and spectroscopy. Given the large number of ongoing photometric and spectroscopic SN surveys, the number of SN events suitable for our method is likely to grow. This will boost our ability to improve and increase the diversity and size of core collapse SN libraries, of key importance to future SN Ia cosmological studies.

As a demonstration of the spectral library, we used the core collapse templates to simulate core collapse contamination in photometrically-selected SN Ia samples. We tested different core collapse modelling approaches, exploiting the luminosity functions and host extinction distributions currently available in the literature. Our analysis suggests that predictions from simulations are sensitive to how global core collapse SN properties (such as luminosity, colour and extinction distributions) are modelled. This suggests that, along with these improved libraries of core collapse SN templates, it is critical to improve our knowledge of the global populations of core collapse SN events.

Acknowledgements

We thank Rick Kessler and Dan Scolnic for their useful comments and feedback. This work was supported by the Science & Technology Facilities Council (STFC) (grant number ST/P006760/1) through the DISCnet Centre for Doctoral Training. MS acknowledges funding from the STFC (grant number ST/N002504/1) and from EU/FP7-ERC grant no. [615929]. This work made use of WISeREP (https://wiserep.weizmann.ac.il) and the Open Supernova Catalog (https://sne.space). This research has made use of the CfA Supernova Archive, which is funded in part by the National Science Foundation through grant AST 0907903. This research made use of Astropy,666http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013; Price-Whelan et al., 2018).

Appendix A Host galaxy extinction estimates and distributions

A.1 Host galaxy extinction corrections

Estimating the extinction on SNe due to dust in the SN host galaxy is challenging. While extinction due to dust in the Milky Way is straight forward to estimate using published galactic dust maps, measuring extinction for SNe located in other galaxies requires more indirect techniques and often relies on empirically-derived relations. Despite the large uncertainties – and indeed possible biases – associated with such host extinction corrections, we still choose to provide host extinction-corrected templates as an option for the end-user. Several SNe included in our library are highly-reddened events, and using these templates in simulations without extinction corrections may not be useful (or at least potentially misleading). In this Appendix, we describe our approach and implementation to extinction corrections.

For most (63 out of 67) of the SNe considered in this paper, we have found a published estimate of the host galaxy reddening, and we therefore apply extinction corrections based on these estimates to provide templates that are corrected for dust at least to first order. These estimates are reported in Table 2 as the colour excess due to host galaxy dust, $E(B-V)_{\mathrm{host}}$.

Most of these estimates of $E(B-V)_{\mathrm{host}}$ have been estimated using the equivalent width of the Na i D absorption lines, and applying the empirically-derived relation presented in Turatto et al. (2003). However, this method presents a large scatter and generally requires high resolution spectroscopy (Barbon et al., 1990; Poznanski et al., 2012; Stritzinger et al., 2018b). An alternative is to use the colour evolution of SNe by comparing the colour evolution of a SN to either intrinsic colour templates (Drout et al., 2011; Taddia et al., 2015), or to samples of SNe thought to suffer minimal extinction (Stritzinger et al., 2018b). This approach is popular in studies of stripped envelope SNe that preferentially occur in regions of high star formation where extinction can be higher and thus more important to correct for. When extinction estimates from this approach are available, we use these in preference to those based on Na i D absorption. For SNe that present negligible reddening ($<0.01$ mag), we set $E(B-V)_{\mathrm{host}}$ to zero. For SNe that lack any information about host reddening (four out of 67 SNe), we use the average $E(B-V)_{\mathrm{host}}$ estimated by Prentice et al. (2016).

In Fig. 11 and Fig. 12 we show the synthetic $B-V$, $V-R$ and $uvw1-V$ colour evolution measured from our templates both before and after our host extinction corrections. For both stripped envelope SNe and hydrogen-rich SNe, the scatter in colour at fixed phase is reduced after extinction corrections are applied, with the mean $B-V$ colour decreased by $0.1$-$0.15$ mag.

We compare these findings with previous analyses. Drout et al. (2011) estimated the scatter in $(V-R)$ at $10$ days after $V$-band maximum ($\sigma_{V-R}$) for a sample of 10 well-observed stripped envelope SNe (8 are included in our library), and find $\sigma_{V-R}=0.06$ mag. We measure $\sigma_{V-R}=0.12$ mag for our 37 stripped envelope SNe (reducing to $0.07$ when only the 8 SNe included in Drout et al. (2011) analysis are considered). If we compare the scatter of $(B-V)$ at $10$ days after $V$-band maximum measured from our stripped envelope SN templates with the sample of minimally-reddened stripped envelope SNe in Stritzinger et al. (2018b) (Table 2), the scatter in our templates is around a factor of two or more larger. These discrepancies are likely due to limitations in the host extinction corrections, but also to the larger diversity included in our templates.

A.1.1 Consequences of incorrect host extinction corrections in the near-UV

Fig. 11 and Fig. 12 suggest that extinction corrections reduce the scatter in the colour of SN templates within the same class. However, some important caveats related to dust extinction corrections need to be highlighted.

Stritzinger et al. (2018b) suggest that SNe Ic preferentially occur in dust regions characterised by a larger value of $R_{V}^{\mathrm{host}}$. Trends like this are more important for our study than the accuracy of an individual host extinction estimate, particularly in the colour of the templates and particularly in the near-UV. In Fig. 13, we show how the extinction changes in the near-UV relative to the optical for different $E(B-V)$ and $R_{V}$. Assuming that the correct value of $E(B-V)_{\mathrm{host}}$ is known, extinction corrections that assume an smaller values of $R_{V}^{\mathrm{host}}$ than the sample average will under-predict the amount of extinction in the near-UV. One consequence of this, from the perspective of using these templates for simulating contamination in SNe Ia cosmological analyses, would be an overestimate/underestimate of the SN Ic contamination at higher redshift where the rest-frame near-UV is redshifted into the optical. For this reason, we tested colour corrections applying the set of $R_{V}^{\mathrm{host}}$ as suggested in Stritzinger et al. (2018b) for different types of stripped envelope SNe: $R_{V}^{\mathrm{host}}=1.1$ for SNe IIb, $R_{V}^{\mathrm{host}}=2.6$ for SNe Ib, and $R_{V}^{\mathrm{host}}=4.3$ for SNe Ic. However we do not observe any reduction in the scatter in the colour evolution of the templates.

A.2 Distribution of host galaxy extinction compared with the literature

In Table 2, we summarise the values of Milky Way and host galaxy reddening assumed in this work. We compare the distribution of host galaxy extinction observed in our sample with the results presented in Prentice et al. (2016) and Hatano, Branch & Deaton (1998). Both studies estimate the distribution of host galaxy extinction for core collapse SNe, the first from a data-driven perspective, the second using simulations.

Prentice et al. (2016) analysed a sample of 110 core collapse SNe with well-constrained host galaxy extinction estimates, and measured the mean and median reddening observed for SNe Ib, Ic/Ic-BL, and II (including SNe IIb and IIn) separately. On average, SNe Ib and Ic/Ic-BL are (respectively) three and two times more extincted than hydrogen-rich core collapse SNe, suggesting differences in the properties of the environments where these events occur.

A different approach was taken by Hatano et al. (1998), who used Monte Carlo techniques to estimate the extinction distribution expected given a set of randomly inclined spiral galaxies, applying some basic modelling of how dust and the progenitors of SNe are spatially distributed in galaxies. They neglected any differences in the spatial distribution of different core collapse SN sub-types. The distribution of $V$-band extinction found by Hatano et al. (1998)777Hatano et al. (1998) presents the distribution of $B$-band extinction, $f(A_{B})dA_{B}$. This can be converted to $f(A_{V})dA_{V}$ assuming the relation $A_{V}=(3.1/4.1)\cdot A_{B}$ can be approximated with the expression:

[TABLE]

This approximate distribution can be implemented in snana, and used to simulate the host galaxy extinction of core collapse SNe.

These two studies are complementary: the first reveals the sub-type dependency of the host extinction distribution, while the second provides an unbiased estimate of the host extinction distribution for core collapse SNe of all sub-types. In Fig. 14 we show the distribution of $E(B-V)$ measured for our set of templates and compare it with the results presented in Prentice et al. (2016) and Hatano et al. (1998). The upper left panel of Fig. 14 shows that the distribution observed in the data is in overall agreement with the predictions from Hatano et al. (1998).

We note the fact that SNe Ib and Ic/Ic-BL present larger tails of highly-reddened objects compared to SNe II, which suggests that different host extinction distributions should perhaps be used for different core collapse sub-types. However, the median extinction measured from our sample of SNe Ib and Ic/Ic-BL are respectively a factor of two and a factor of four lower than that in Prentice et al. (2016). This suggests that different biases affect the samples considered, with our sample likely excluding some highly-extincted objects.