A statistical approach to identify superluminous supernovae and probe their diversity

We begin by constructing our SLSN I sample from published events in the literature. We require a well-observed sample, with at least six epochs of photometric data between phases of $-15$ d to $+30$ d, of which at least one must lie between $-15$ d and $0$ d, sampling a synthetic box filter at rest-frame 4000 Å. This synthetic filter was introduced in Inserra & Smartt (2014), together with a second box filter at 5200 Å, which are designed to sample two regions of SLSN I spectra that are dominated by continuum, with few absorption features. In this paper, all phases are given in rest-frame days relative to peak brightness in our synthetic 400 nm filter.

We further select our SLSN I sample to only contain events with a single main peak in the light curve, as otherwise there can be ambiguity in identifying the main peak and measuring phases. We therefore exclude those events with a ‘secondary peak’ (e.g. iPTF13dcc, iPTF15esb; Vreeswijk et al., 2017; Yan et al., 2017b, around 5% of the literature SLSNe I), where a ‘secondary peak’ is defined as a peak, before or after the brightest (main) peak, showing an absolute magnitude difference of $\lesssim 1$ mag with respect to the main peak. We also remove events with H$\alpha$ in their spectra (e.g., iPTF13ehe; Yan et al., 2015), another 5% of SLSNe I in the literature. Note we do not exclude a priori SLSNe I with early-time ‘bumps’ (see Table 1).

Around 50% of the published SLSN I events in the literature pass our requirements, and these can be found in Table 2, while those rejected are reported in Table 1 together with the criterion by which they are removed from the sample.

We $k$-correct all published photometry for each object to our two synthetic filters (400 nm and 520 nm), calculated with the snake¹¹1https://github.com/cinserra/S3 software package (Inserra et al., 2016b), which also estimates the uncertainties on the $k$-corrections. The synthetic filters, together with the standard Bessell $B$/$V$ filters and Sloan Digital Sky Survey (SDSS) $u$/$g$ filters, are shown in Figure 1. Applying $k$-corrections from arbitrary observed filters over $0.1<z<4.0$ to $B$ and $V$ instead of the two box filters would result in a difference of $-0.03$ mag at peak epoch (both $400-B$ and $520-V$), and $0.01$ mag ($400-B$) and $0.05$ mag ($520-V$) around 30 days after rest-frame maximum. When observed spectra for a specific SLSN are not available, we use an average SLSN I time-series spectral energy distribution (SED), based on the methodology of Prajs et al. (2017). We correct all our observed photometry for Milky Way extinction prior to $k$-correction using the prescriptions of Schlafly & Finkbeiner (2011), but make no corrections for extinction in the SN host galaxies, which is believed to be small (e.g. Nicholl et al., 2015b; Leloudas et al., 2015b). Finally, we convert the rest-frame apparent magnitudes into absolute magnitudes using a flat $\Lambda$CDM cosmology, with $H_{0}=72$ km s^-1 Mpc^-1, $\Omega_{\mathrm{matter}}=0.27$, and $\Omega_{\Lambda}=0.73$.

To estimate the SLSN I brightness around peak epoch ($-15\leq\mathrm{phase}\leq 30$), where literature SLSNe I have the most coverage, we investigate several techniques to fit the available data set, including polynomial fitting (as in Inserra & Smartt, 2014), and interpolation using Gaussian processes (GPs) regression (Bishop, 2006; Rasmussen & Williams, 2006). GPs are already successfully used in several areas of astronomy (e.g. Mahabal et al., 2008; Way et al., 2009; Gibson et al., 2012) and in the context of supernovae (e.g. Kim et al., 2013; Scalzo et al., 2014; de Jaeger et al., 2017). In supernova analyses, GPs can be used for Bayesian regression and mean function fitting with a non-parametric approach, e.g. broad-band light-curves, bolometric light-curve, temperature and radius evolution, as well as line profiles in spectra.

A GP assumes that our variable $y$ is randomly drawn from a Gaussian distribution with a certain mean and covariance ($y\sim f(\mu,\sigma^{2})$), and then considers $N$ such variables drawn from a multivariate Gaussian distribution computing their joint probability density. This is $Y\sim f(m,K)$, where $Y=(y_{1},...,y_{N})^{T}$, $m=(\mu_{1},...,\mu_{N})^{T}$ is the mean vector (transposed) and $K$ is the covariance matrix, called the ‘kernel’, having a $N\times N$ dimension of the form

In the case of uncorrelated variables $cov(y_{i},y_{j})=0$, this becomes a diagonal matrix. This approach provides for a probability distribution over functions and it allows us to compute a confidence region for the underlying model of our variable. A GP is hence specified by its mean function and kernel.

To reach convergence of the distribution, the kernel hyper-parameters are optimized using the maximum likelihood method. We test several kernels to find the most suitable covariance function for the objects in our SLSN sample. The kernels we consider are an exponential sine squared kernel (suited for periodic functions), a linear kernel, a polynomial kernel and a squared exponential kernel. As a basic metric for quality of fit we compute the $\chi^{2}$, comparing the mean of the GP posterior distribution to our data. We find that the ‘Matern-3/2’ kernel gives the best fit. The kernel can be written in terms of radius, $r=|a_{i}-a_{j}|$, as:

where $c$ and $t$ are the best fit hyper-parameters.

We use the Python package george (Ambikasaran et al., 2014) to perform our GP regression. Figure 2 (top) shows a comparison between the Matern-3/2 kernel and a third-order polynomial fit. The GP fit is a better representation of the data with respect to the polynomial, meaning a lower $\chi^{2}$. Figure 2 also shows examples of the GP fits to 400 nm data, while all the fits are reported in Appendix A. As expected, when the data are sparse the fit uncertainties are larger (bottom panel). A key advantage of the GP fitting is that the fit uncertainties, as a function of phase, are naturally produced by the GP fitting. Accurate uncertainties will be important for the analysis in this paper. We further refer to Ivezic et al. (2014) for a more in-depth analysis of advantages and drawbacks of GPs in astronomy.

We fit our entire SLSN sample from -15 d (in the 400 nm band) to +55 d. Whenever the data or GP uncertainties on the magnitude are smaller than those reported by the survey that discovered the SLSN, we replace them with the typical survey photometric uncertainties at the redshift of the SN (e.g., for the 400 nm band, PS1 averaged uncertainties at $z<0.25$ and $0.25<z<0.60$ are 0.02 and 0.06 mag, while the DES uncertainties at $z<0.60$ and $0.60<z<0.90$ are 0.04 and 0.05 mag)²²2Derived from SLSN PS1 and DES papers listed in Tables 2 & 1 and from DES private communications.. The resulting fit magnitudes are reported in Table 2, together with their uncertainties.

Our final measurements concern the SLSN spectra, and in particular the estimation of the photospheric velocities. In core collapse SNe, Sc II $\lambda$6246, and subsequently Fe II $\lambda$5169, are the best available proxies to trace the photospheric evolution due to their small optical depth (Branch et al., 2002). In SLSN I spectra, Sc II is not visible, but Fe II has been measured using several different approaches (e.g. Inserra et al., 2013; Nicholl et al., 2015b; Liu et al., 2016).

We measure the line velocities in all spectra from +10 to +30 d. Before +10 d, the ionic component is weak (Inserra et al., 2013) and contaminated by Fe III (Liu et al., 2016). When spectra around $10\pm 2$ d and $30\pm 2$ d are not available for a given SLSN, we estimate the velocity from a least-squares fit of the measurements from nearby epochs and account for an additional uncertainty in the estimate using $\sigma_{\rm final}=\Delta t/\sigma_{\rm measure}$, where $\Delta t$ is the phase difference (in seconds) between the estimated and measured epochs. Only 12 out of the 19 SLSNe have spectra covering the wavelength region and the time-frame of interest (see Table 2).

We measure the velocity from fits to the absorption minima. We experiment with three different profiles for the fits (Gaussian, skewed Gaussian, and Voight), finding an overall agreement among the three profiles. We repeat the measurements several times for each feature, changing the continuum levels to better estimate the uncertainties. We then use the mean of the measurements as the final value, and the standard deviations as the uncertainty estimate; the values are tabulated in Table 2. This approach has been widely used in measuring the line velocities of SLSNe I, with consistent results (e.g. Pastorello et al., 2010; Chomiuk et al., 2011; Inserra et al., 2013; Smith et al., 2017).

We have also cross-checked our velocity measurements using a GP approach (Section 2.2), fitting the wavelength region from 4800 Å to 5200 Å, and finding the local minimum. We use the same Matern-3/2 kernel as with the light curve fitting. As would be expected, given the greater flexibility, the GP fits usually return a $\chi^{2}$ comparable to or lower than the standard fitting procedure, with the advantage of improving fitting/deblending multi-component profiles without a biased prior knowledge of the feature types and numbers (e.g., absorptions/emissions with P-Cygni profiles, Lorentzian or Voight wings, etc.). However, to properly evaluate the uncertainties we need a kernel function based on the uncertainties in the flux of the spectra. This information is missing for $\sim$40% of our data set, and hence we use the profile fitting to allow for consistency in the approach.

Our measurements and line evolutions are broadly in agreement with those of Liu et al. (2016), those reported in the papers listed in Table 2, and the photospheric velocities reported by the modeling of Mazzali et al. (2016). The only noticeable difference is in the velocity of Gaia16apd, where we found a decrease of $\sim$1000 km s^-1 over the phase range analyzed. This is due to the presence of galaxy lines that make the fit more complicated and biased by the choice of the number of components to analyze. The Fe II $\lambda$5169 Å velocity measurement is reported as that at +10 d, or $v_{\mathrm{10}}$ (km s^-1), while the velocity evolution over the phase range from 10 d to 30 d post-peak is $\dot{v}=-\Delta v/\Delta t$ (km s^-1 day^-1), in a similar fashion to that used in SNe Ia (Benetti et al., 2005).

The parameter space of Figure 3 may in principle be used to help physically understand the explosion mechanisms of these transients. Relationships between the change in luminosity in one band (panel A) and the color evolution (panel B), and the broad band behavior of a SN at a given epoch (panels C & D) are broad reflections of the physical properties of the SN ejecta (i.e. diffusion time, opacity and temperature). The correlation shown within panel A is likely a reflection of the diffusion time of the ejecta – similar to that seen within SNe Ia (Phillips, 1993). However, as we do not consider the light curves of SLSNe I to be radioactively driven (e.g. Nicholl et al., 2013; Inserra et al., 2017a), for SLSNe I this correlation is unlikely to be solely related to the mass of the ejecta produced.

On the other hand, the tight relation presented in panel D of the 4OPS (see Figure 3 and Table 3) between the color observed at peak and at +30 d suggests that these two are correlated by one physical parameter only, which could be the temperature or the radius.

A wide range of possibilities have been postulated to explain SLSN I luminosities, such as the rapid spin-down of a magnetar (e.g. Kasen & Bildsten, 2010; Woosley, 2010; Dessart et al., 2012), the interaction between the SN ejecta and the surrounding CSM previously ejected from the massive central star (e.g. Chatzopoulos et al., 2012; Woosley, 2017), and a pair instability explosion (e.g. Kozyreva et al., 2017). For all three models there are multiple parameters at play in the production of the overall luminosity and color evolution of the transient. As such, the linking of an observed behavior to a dominant physical parameter becomes complex. For example, a magnetar magnetic field strength, spin-period and explosion ejecta mass are all factors in the luminosity evolution (Kasen & Bildsten, 2010), whilst within the interaction model, the mass of the ejecta, its density profile and distribution coupled with the mass, distance and volume of the CSM shell must be considered (e.g. Chatzopoulos et al., 2013; Woosley, 2017).

At present there are no model predictions that aptly describe the broad-band behavior shown in Figure 3. This could be due to the fact that the diffusion time not only depends upon the ejected mass, but also on the ejecta velocity and its opacity to optical-wavelength photons. Opacities, in particular, are determined by the temperature and composition of the ejecta and therefore may vary with time during the SLSN I evolution (Mazzali et al., 2016). Thus, explaining the relation observed within any of panels A, B and C of the 4OPS with any of the above suggested scenarios is not trivial. Nonetheless, the presence of the relations hints that a pure radiative transfer in the SN ejecta should be at play and any model that aims to explain such SNe should take these observational properties into account.

The primary purpose of our work is to define a main population of SLSNe I. Within the context of a magnetar powered event, favored by several observational studies (e.g. Chen et al., 2015; Nicholl et al., 2015b; Smith et al., 2016; Inserra et al., 2016a), the behavior observed could be explained by the magnetar energy injection always occurring at a certain time (e.g., shortly after the explosion). Such a scenario would allow the diffusion time of the ejecta to be comparable to the time needed for the SN to reach peak luminosity (i.e., panel A). An injection this early would also provide the rotational energy needed to overwhelm the initial thermal energy of the SN explosion, and hence provide the energy source which drives the main peak of the light curve (Kasen & Bildsten, 2010). Such a population of engine driven SLSNe I would be composed of brighter objects that are overall bluer and more slowly evolving than dimmer events. Moreover, redder objects in this main population would evolve faster in both luminosity and color as inferred from Figure 3, panel D, and previously shown by Inserra & Smartt (2014).

Objects outside the main population of SLSNe I, but with a luminosity evolution that could be explained by an inner engine, have already been found (e.g. Greiner et al., 2015; Kann et al., 2016, and the Dark Energy Survey collaboration, private communications), but their spectrophotometric evolution is different to SLSNe I, and this is reflected in Figure 5. These objects could belong to a similar engine-driven transient family, only here the injected energy and/or the timescale over which it is injected would be somewhat different.

The analysis of Section 4 returns two subclasses, which are outlined in terms of spectrophotometric evolution during the first 30 days from peak, as well as a distinct photospheric velocity behavior. In the context of the magnetar scenario, the almost flat velocity evolution exhibited by the slow subclass, and their overall slower velocity, suggests that the photosphere reaches the internal shell created by the magnetar bubble (see eq. 7 in Kasen & Bildsten, 2010) earlier than in the fast subclass. Fast SLSNe I with a high photospheric velocity also have a larger velocity gradient. This could be related to additional energy deposited by the magnetar into the ejecta (Dessart et al., 2012). Such energy is a function of the spin period (see eq. 1 in Kasen & Bildsten, 2010), and faster rotation would imply more energy and hence faster ejecta.

The almost frozen spectral evolution exhibited after peak by slow SLSNe I (Nicholl et al., 2016; Inserra et al., 2017a) supports our findings and the idea that the photosphere reaches the inner shell earlier than in the fast events. In addition, the slow events show forbidden lines earlier, suggesting they become optically thin earlier. This could be explained by a high amount of oxygen ($\sim 10$M_⊙, Jerkstrand et al., 2017), whose recombination could hasten the process.

Generally, differences in properties in SN subclasses, and hence between the slow and the fast, could be due to different degrees of mixing or geometry (Leloudas et al., 2015a; Inserra et al., 2016a; Leloudas et al., 2017), which is true regardless of the source of additional luminosity. The photospheric velocity depends on the the optical depth, for which the heavy elements with higher line opacities are the prime contributor. A more efficient mixing of heavy elements in the outer ejecta might result in an initial higher photospheric velocity, whereas a less efficient one could lead to slow velocity and constant temperature. The gradient evolution may also be explained in terms of the ejecta density structure or of the photosphere moving to non-mixed layers of the ejecta.

Figure 3 can be used to define a homogeneous population of events, with 90% of previously classified SLSNe I meeting our definition, for further study in a cosmological context. This can be used for current (e.g,. DES, Pan-STARRS) and new generation (e.g., LSST, Euclid and WFIRST) surveys to identify/classify (also in real time) SLSNe I. This would allow identification even without a spectroscopic evolution, important given relatively limited spectroscopic resources. Moreover, with only a spectrum at +10 days, the identification can be confirmed as well as distinguishing between the fast and slow subgroups. This, together with Figure 3, strengthens further the possibility to use them as standardizable candles at high-redshift. The next logical step is to discern the two subclasses only with photometry and/or to move this analysis to shorter wavelengths (e.g., the rest-frame ultraviolet) allowing higher redshifts to be studied.