The Hydrogen-poor Superluminous Supernovae from the Zwicky Transient Facility Phase I Survey. I. Light Curves and Measurements

Here, we briefly describe the photometric selection system of our SLSN-I candidates. Daily ZTF alerts are ingested into a dynamic science portal called the GROWTH Marshal (Kasliwal et al. 2019). A filter is implemented within the Marshal to select SLSN candidates. The candidates passing the filter are not automatically saved; instead, each week, a human scanner has to visually examine the LC of each candidate selected by the filter and make a decision as to whether it is worth being saved. The subsequent spectral follow up is based on the human-saved candidates. This filter adopts several cutoff conditions, including: (1) it is not moving—the same alert is detected in two consecutive epochs with Δt > 0.02 day; (2) it is not a star, based on the Sloan Digital Sky Survey (SDSS) star–galaxy scores (Tachibana & Miller 2018); (3) the exclusion of bogus, alerts based on the scores constructed by Duev et al. (2019); (4) it is not in the galactic plane, with galactic latitude ∣b∣ > 7°; and (5) variability has not been detected at this location more than a year prior to the alert.

In addition, we assign numerical scores to several properties, including: (1) slowly rising events; (2) faint blue hosts; (3) the spatial location of the transient relative to the center of the host (against nuclear transients); (4) the brightness of the candidate relative to the host brightness (against faint transients); and (5) the time interval between the first detection and the last one (against very long-lived transients, like active galactic nuclei). This candidate filter gives high scores (thus preferences) to slowly rising events with faint blue hosts. For example, a candidate rising at least 20–25 days, with a faint host galaxy, will have a high score, pass the filter, and also be saved by the human scanner.

We note that the ZTF collaboration has multiple transient groups that also perform daily alert stream scanning, and most transient candidates can be saved by multiple groups—for example, the ZTF Bright Transient Survey (BTS; Fremling et al. 2020), the ZTF Census of the Local Universe (De et al. 2020), the fast transient group, the stripped envelope SN group, and the ZTF nuclear transient group. Therefore, the classification of a specific SLSN-I candidate can also be drawn from other classification efforts, most noticeably the BTS. A few classifications are taken from the Transient Name Server reported by external groups. But almost all of these externally classified sources were also identified as candidates by our filter. The filtering method will be described in detail in the forthcoming work.

With the above selection criteria, on the order of 50 candidates per week are saved, before going through another round of vetting when the spectral classification observations are planned. Any candidates brighter than 18.5 mag are classified by the ZTF BTS. The classification efficiency at ≤18.5 mag is very high, close to 95% (Fremling et al. 2020; Perley et al. 2020). We anticipate that our catalog is almost entirely complete for SLSNe within the ZTF footprint peaking at magnitudes brighter than 18.5 mag, and highly complete to 19.0 mag. For magnitudes fainter than 19.0 mag (44 of 78 events in total), there may be biases, primarily related to the nature of the host and/or the rising phase of the LC, which will be examined in future work. The detailed structure of the LC, e.g., the presence of bumps, was not a crucial factor for selection or triggering follow up.

Our efforts of spectroscopic classification were primarily focused on SLSN candidates fainter than 18.5 mag, using the Spectral Energy Distribution Machine (SEDM; Blagorodnova et al. 2018) and the Double Beam Spectrograph (DBSP; Oke & Gunn 1983) mounted on the Palomar 60 inch (P60) and 200 inch (P200) telescopes, respectively. Additional facilities include the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on the Keck I telescope, the Alhambra Faint Object Spectrograph and Camera (ALFOSC) on the 2.56 m Nordic Optical Telescope (NOT), the SPectrograph for the Rapid Acquisition of Transients (SPRAT) on the 2 m Liverpool Telescope (LT), and the Intermediate-dispersion Spectrograph and Imaging System (ISIS) on the 4.2 m William Herschel Telescope (WHT). The basic information about the classification spectra is listed in Table A2. The spectral reduction is performed using various standard reduction pipelines. This includes the SEDM automated pipeline (Rigault et al. 2019), the pyraf-dbsp package (Bellm & Sesar 2016) and DBSP_DRP (Roberson et al. 2022) pipelines for the DBSP data, and the LPipe package (Perley et al. 2019) for the LRIS data.

Finally, we note that the SLSN project has at least 0.5–1 night of DBSP time on P200 per month for spectral classification (PI: Yan). More quantitative discussions on the completeness of the spectral classification will be included in Paper III.

Every event in our sample has at least one spectrum, and most have multi-epoch spectra. The hallmark spectral features for SLSNe-I are the five O ii absorption features in the wavelength range of 3737−4650 Å in prepeak and/or near-peak optical spectra. These were identified and discussed in Quimby et al. (2011, 2018). We utilize the large SLSN-I spectral template library assembled by Quimby et al. (2018), and update the library by adding the missing phase information to some of the templates. Our classification relies on matching with the spectral templates using SNID (Blondin & Tonry 2007) and superfit.

To determine the best-matched spectral templates, we run both superfit and SNID on the smoothed spectra, with host galaxy emission lines removed and redshifts fixed for most sources. If the spectrum of an SLSN-I has a good match with that of an SN Ic, but its g-band peak luminosity is higher than −20.0 mag, we classify this candidate as an SLSN-I, since many SLSNe-I develop spectra similar to those of SNe Ic after the peak (Pastorello et al. 2010; Quimby et al. 2018; Gal-Yam 2019).

In the Appendix, Figure A1 presents the best-fit spectral templates for each event in our sample, with the event name, phase, and template information labeled after the spectra. The phases are measured relative to the rest-frame g-band peak in this paper. We note that the phase differences between the observed spectra of our SLSNe-I and the templates in the library are not zero, but generally within 30 days (rest frame). This is not surprising, since SLSNe-I can have similar photospheric spectra, but different LC evolution timescales (see Kangas et al. 2017; Quimby et al. 2018).

The above classification procedure can sometimes give ambiguous results for a small number of events—i.e., two events for our sample. In these cases, we rely on additional LC information, such as the rise time and peak luminosity, to break the degeneracy. For example, ZTF19aacxrab (SN 2019J) and ZTF19aaqrime (SN 2019kwt) are almost equally well matched with the spectral templates of SN Ia and SLSN-I. In the case of SN 2019J, although the overall spectral features broadly match with those of SNe Ia, its spectrum lacks S ii λλ 5433, 5606 commonly seen in SNe Ia. In addition, because of their slow-rising LCs and high peak luminosities (i.e., SN 2019kwt has M_g ∼ −22.8 mag), we adopt the SLSN-I classifications corresponding to the template spectra of SN 2007bi and PTF09cnd, respectively.

In summary, the 78 sources listed in Table A1 can be classified as SLSNe-I according to the features of their spectra and LCs. The classification spectra are made available to the public as part of the electronic data at the Journal website, and they will be uploaded to the Weizmann Interactive Supernova Data Repository (Yaron & Gal-Yam 2012). ²⁷

The bulk of the photometric data comes from ZTF, including the data from the public survey with a 3 day cadence and the ZTF partnership and Caltech surveys with a faster cadence (≤2 days) over smaller areas (Bellm et al. 2019b).

The IPAC ZTF pipeline produces reference-subtracted images using the ZOGY algorithm (Zackay et al. 2016) and aperture photometry for all transients detected at ≥5σ. However, there are two issues with the LCs produced by this pipeline. One is that the reference images built in 2018 may contain signals from transients that exploded during the same period, for which the photometric offset needs to be corrected. The second issue is that the upper limits prior to the first detection are from aperture photometry and based on the image noise measured over the entire quadrant. This can significantly underestimate the transient signals.

To fix these two issues, we perform forced point-spread function (PSF) photometry using the software provided by IPAC. ²⁸ With a code from Yao et al. (2019), we refine the astrometric position of each event by using only the images around the peak phase. The very early- and late-time forced photometry without transient signals allows us to compute baseline offsets to the LCs. In addition, we reject bad-quality data if: (1) the image processing and instrumental calibration fail to meet predefined quality criteria; (2) the robust estimate of the 1σ value of the spatial noise per pixel in the image is over 25; and (3) the seeing is larger than 5'', similar to what was used in Yao et al. (2019). The photometry of the same transient observed with different CCD quadrants can have a systematic offset. This problem is fixed by our photometric reprocessing as well. For the final photometric collections, a detection is defined as 4σ (3σ for Swift data) above the background, and an upper limit is computed at 3σ.

The ZTF astrometric and photometric systems are calibrated using the Gaia Data Release 1 data and the PS1 catalogs, respectively (for details, see Masci et al. 2019). The output magnitudes are in the AB system. The airmass and color term corrections are included when the photometry is calibrated to the PS1 system. The color information for color corrections (g − r for the g, r bands and r − i for the i band) is obtained from the g, r, i photometry taken at the same epoch, or neighboring observations.

3.2.1. P60, LT, and P200 Photometry

3.2.2. Swift Data

4.1.1. Extinction Corrections

The galactic reddening E(B − V) is taken from the Schlafly & Finkbeiner (2011) dust map, using the NASA/IPAC InfraRed Science Archive database. The extinction corrections at different wavelengths are computed using the empirical dust extinction laws of Fitzpatrick & Massa (2007), with R_V = 3.1.

Previous studies have shown that SLSN-I hosts are mostly low-mass metal-poor dwarf galaxies (Lunnan et al. 2014, 2015; Leloudas et al. 2015; Angus et al. 2016; Perley et al. 2016; Chen et al. 2017a; Schulze et al. 2018). This is the case for most of the events in our sample. Because their LCs do not show particularly red colors at prepeak phases, we do not make any host galaxy reddening corrections to these events. However, a small fraction of the SLSN-I hosts are massive (10⁹⁻¹⁰ M_⊙; Perley et al. 2016; Chen et al. 2017b), and the host galaxy reddening in such cases can be large. We find that seven events in our sample (Table 2) have non-negligible host galaxy reddening. They either reside in bright hosts (M_r ≲ −18.5 mag) or have redder (g − r) colors at peak than average (see Figure 8), which are presumably due to host galaxy reddening.

Table 2. Host Galaxy Reddening

Name	$E{\left(B-V\right)}_{\mathrm{host}}$ (mag)	Method
SN 2018don	0.4 a	spec.temp
SN 2018kyt	0.11	line ratio
SN 2019kwt	0.22	line ratio
SN 2019aamx	0.07	spec.temp
SN 2019aamu	0.23	spec.temp
SN 2019stc	0.18	line ratio
SN 2020xkv	0.24	spec.temp

Note.

^a Estimated in Lunnan et al. (2020).

Download table as:  ASCIITypeset image

The host galaxy reddening is difficult to measure accurately, but it can be roughly estimated using two different methods. The first one uses the Balmer line ratio measured from the host galaxy spectra. Without dust attenuation, the predicted line ratio between Hα and Hβ remains constant over a wide range of temperatures and electron densities in the nebular regions (e.g., Hα/Hβ ≈ 2.86; Osterbrock & Ferland 2006). The attenuation caused by dust scattering is wavelength-dependent and stronger at shorter wavelengths. This method is affected by both stellar and SN absorption features. We remove the stellar Balmer absorption using FIREFLY (Wilkinson et al. 2017), before we measure the line flux, but we ignore the influence of SN features. This method is applied to the host spectra of ZTF18acyxnyw (SN 2018kyt) and ZTF19acbonaa (SN 2019stc). Note that the SN location may not coincide with the H ii region responsible for the Balmer lines, so the real extinction at the SN site could be different from the inferred value. But such an estimation is at least an indicator of whether the galaxy has significant dust.

The second method is to infer the host galaxy reddening by matching the observed spectrum (including the spectral continuum slope and absorption features) with an SLSN-I template with negligible extinction and at a similar phase. For this method, the assumption is that the observed red spectral color is due to the host galaxy reddening. The results from this method are very uncertain, but offer crude extinction estimates for those events without host emission lines. This method is applied to ZTF18aajqcue (SN 2018don), SN 2019aamx, and ZTF19acvxquk (SN 2019aamu), whose spectra lack distinct Balmer emission lines from host galaxies. For example, to match the near-peak spectrum of SN 2018don, the t ∼ + 54 day template spectrum of SN 2007bi is required to be reddened by E(B − V) ∼ 0.4 mag (Lunnan et al. 2020). Similarly, the +9 day spectrum of SN 2019aamu matches well with the −7 day template spectrum of PTF12gty, after considering a host galaxy reddening of E(B − V) = 0.23 mag. And the +13 day spectrum of SN 2019aamx matches well with the +22 day template spectrum of PTF09cwl with E(B − V) =0.07 mag (shown in Figure A1). Note that both PTF12gty and PTF09cwl have a faint dwarf host, so their host galaxy reddenings are assumed to be negligible (Perley et al. 2016; Quimby et al. 2018).

For SN 2019kwt and ZTF20abzaacf (SN 2020xkv), we applied both of the above methods, since their spectra contain strong SN signals and we are not able to remove the stellar absorption or SN features. The derived E(B − V) for SN 2019kwt is 0.22 mag from the host galaxy spectrum and 0.18 mag from template matching, which appear to be consistent. In the case of SN 2020xkv, the E(B − V) derived from the spectral matching ranges from 0.22–0.26 mag, in contrast to the lower value (0.17 mag) inferred from the Balmer decrement. Its spectrum has a very low signal-to-noise ratio at Hβ, and with both being highly uncertain measurements, we adopt the average value of 0.24 mag from the template matching method.

4.1.2. The K-corrections

The absolute g-band peak magnitude is one of the key parameters of SLSNe-I. Proper estimates of this value require proper K-corrections to derive reliable rest-frame values, either in the g band for sources at z ≤ 0.17 or in the r band for sources at z > 0.17. This switch is chosen because the observed r band LC at z > 0.17 is closer to the rest-frame g band in wavelength.

Most events in our sample have at least one spectrum around the LC peak (i.e., ≲20 days in the rest frame), which allow us to compute the K-corrections using the formula described in Hogg et al. (2002). For those without spectra near the peak phases, we apply a constant K-correction of $-2.5\,\times \mathrm{log}(1+z)$. This approximation is not too far from the spectral estimates, as shown below.

Figure 2 compares the calculated K-corrections based on the observed spectra (blue dots) and an assumed 10⁴ K blackbody spectral energy distribution (SED; green line) with the constant correction (black line). The constant value of $-2.5\times \mathrm{log}(1+z)$ appears to capture most of the corrections within a range of ±0.1 mag (the shaded region) for SLSNe-I at 0.06 < z < 0.67. We list the K-corrections for each event in Table A4.

Zoom In Zoom Out Reset image size

4.1.3. The S-corrections

As the g-, r-, and i-band photometry of the SLSNe-I presented in this work is obtained from four different telescopes, the errors caused by filter differences need to be evaluated. Since most of the data are from the ZTF itself, we choose to correct the LT/P200/P60 photometry to that of the ZTF filters. Following the method of Stritzinger et al. (2002, 2005), Pignata et al. (2008), and Wang et al. (2009), the correction between the different instruments (known as the S-correction) can be computed using

$$\begin{eqnarray}&&{{Sc}}_{\lambda 1}={M}_{\lambda 1}-{m}_{\lambda 1}-{{CT}}_{\lambda 1}\left({m}_{\lambda 1}-{m}_{\lambda 2}\right)-{{ZP}}_{\lambda 1},\end{eqnarray} \tag{ 1 }$$

where M_λ1 is the SN synthetic magnitude computed with the response functions of the ZTF λ1 filter and m_λ1 and m_λ2 are the SN synthetic magnitudes computed with the LT/P200/P60 filters. CT_λ1 is the color term and ZP_λ1 is a zeropoint, which are measured by convolving the ZTF filters with a large sample of spectrophotometric Landolt standard stars from Stritzinger et al. (2005). The response function above includes filter transmission, detector quantum efficiency, and atmospheric transmission.

As S-corrections are spectrum-dependent, and accurate measurements at any given epoch require far more SLSN-I spectra than we have in this work, we use the S-correction values computed from the spectra as an additional 1σ error in photometry, and add them to the existing errors of LT/P200/P60 photometry in quadrature. The additional 1σ errors due to the S-corrections are listed in Table 3.

Table 3. The Estimated Error Caused by the S-correction

Telescope	Filter
	g(mag)	r(mag)	i(mag)
LT	0.018	0.044	0.030
P60	0.042	0.030	0.026
P200	0.023	0.030	0.025

Download table as:  ASCIITypeset image

The LC analysis requires estimates of various parameters—such as the peak magnitude, rise time, and rise rate—which all involve numerical interpolation and fitting. We adopt a machine-learning algorithm, Gaussian Process (GP) regression, which has various kernel functions. The GP interpolation can reduce the influence of outliers and gives robust error estimates. We use a composite kernel, which is the sum of a Matérn kernel with a white noise kernel. We tested Matérn kernels with ν parameters of 3/2 (Matérn 3/2; Inserra et al. 2018a; Angus et al. 2019; Lunnan et al. 2020) as well as ν = 5/2 (Matérn 5/2). When we compute the peak magnitudes and phases, a GP fit with Matérn 3/2 is performed in the flux space. The Python package george (Ambikasaran et al. 2016) and Scikit-learn give comparable results.

To set an approximate luminosity scale, we also plot the absolute magnitude on the right Y-axis of Figures A2–A7. This is calculated by assuming a constant K-correction of $-2.5\times \mathrm{log}(1+z)$. In these figures, the LCs have been corrected for galactic and host galaxy reddening, whenever possible.

To calculate the apparent peak magnitude and phase, we run the GP to interpolate the LCs in flux space. The errors are shown as 1σ uncertainties. The g-band absolute peak magnitudes are derived using M_g = m_g − μ − KC, where μ is the distance modulus, KC is the K-correction, and m_g is the apparent g-band magnitudes, corrected for both the host (in the rest frame) and galactic extinction (in the observed frame). The rest-frame g-band absolute peak magnitudes and the peak dates are tabulated in Table A4. The peaks of some events are not well constrained, due to the poor sampling around the peak. Figure 3 shows the distribution of the g-band absolute peak magnitudes (M_g ) for the ZTF sample, ranging from −19.8 mag to −22.8 mag, with a median and 1σ error of ${M}_{g,\mathrm{med}}=-{21.48}_{-0.61}^{+1.13}$ mag. Both our median value and dispersion are consistent with those of previous samples.

Zoom In Zoom Out Reset image size

The weak bimodal distribution in the peak absolute magnitudes should not be overinterpreted, because this is the raw distribution function, without applying corrections for various selection biases (e.g., the Malmquist bias or the preference of our selection system for slowly rising events). Selection biases could also have a strong impact on the distribution of the timescales (e.g., Figure 5), the peak bolometric luminosities (e.g., Figure 16), and other photometric properties. The correct interpretation of these distributions requires a well-defined selection function for our sample—a detailed simulation of this, however, is well beyond the scope of this paper. Further studies (e.g., Paper III) are required to confirm these distributions.

Traditional rise times, defined as the interval between the explosion date and the peak, are notoriously difficult to measure, especially for distant transients such as SLSNe-I. Although the ZTF detection can go down to 20.5–21.0 mag, at the median redshift of our sample it only probes an absolute magnitude of ∼−20 mag. The upper limits before the first detection cannot usually place strong constraints on the LC evolution in the very early phases, and only about one-third of our sample has high-quality early-time data. In addition, Anderson et al. (2018) show that an SLSN-I (SN 2018bsz) can have a long slowly rising "plateau," before a steeper, faster rise to the peak, which makes it harder to speculate on explosion dates without enough deep early data.

Instead of using the traditional rise time, we define t_rise/decay,x as the time interval between the peak and the flux when it is at a fraction x of the peak value. Here, the x factor can be 10% (Δmag = 2.5), or 1/e (Δmag ≈ 1.09). Time dilation is corrected for all timescale measurements. The derived t_rise,10% and t_{rise/decay,1/e} are listed in Table A4. The poorly constrained timescales are not included in the following analysis.

Figure 4 shows the direct proportionality between t_rise,1/e and t_rise,10% in the rest-frame g band. A linear fit gives t_rise,1/e = 0.80 t_rise,10% − 1.73 days with the 1σ uncertainty as small as ∼2.84 days. The small scatter implies that the LC rise rates in the very early phases do not have significant differences in our sample. The rise times t_rise,10% cover a wide distribution, ranging from 10 to 90 days, with the mean (median) value and a standard deviation of $\overline{{t}_{\mathrm{rise},10 \% }}=41.9(38.3)\pm 17.8$ days.

Zoom In Zoom Out Reset image size

Of the total 56 SLSNe-I with t_rise,10% measurements, the two slowest rising events—SN 2018ibb and ZTF20aapaecd (SN 2020fyq)—stand out in Figure 4, with timescales longer than $\overline{{t}_{\mathrm{rise},10 \% }}+2\sigma =78$ days.

For another slowly evolving event, ZTF20aadzbcf (SN 2020fvm), its t_rise,1/e reaches 91 days, which is suggestive of an unusually long rise time, although the t_rise,10% is not well constrained. SN 2020fvm has two LC peaks and we set the second—and also brighter—one as its main peak. If the first one is set as the main peak, its t_rise,1/e is around 48 days. The fastest event in our sample is SN 2018bgv (Lunnan et al. 2020), which rises to the peak in less than 10 days. Such a fast and luminous event is very rare, and it may be associated with the Fast Blue Optical Transients (Drout et al. 2014; Ho et al. 2021). Compared with previous SLSN-I samples, our sample contains slightly more fast-evolving events. For instance, nine of 56 events (9%–26%, calculated at a confidence level of 95%, based on Gehrels 1986) are found to have t_rise,1/e ≲ 15 days, while only five of 55 events (4%–18%) were reported in the previous samples.

Outside the ±1σ range, we have 10 fast-evolving events with ${t}_{\mathrm{rise},10 \% }\leqslant \overline{{t}_{\mathrm{rise},10 \% }}-1\sigma =24$ days and 10 slow events with ${t}_{\mathrm{rise},10 \% }\geqslant \overline{{t}_{\mathrm{rise},10 \% }}+1\sigma =60$ days. The t_rise,10% shows a continuous distribution, and it cannot be divided into two separate fast and slow subclasses, as indicated by previous studies (Nicholl et al. 2015; De Cia et al. 2018).

Figure 5 shows a comparison of the rise and decay timescales measured at 1/e of the peak flux for 48 events. Both the rise/decay timescales show a distribution centered at ∼25–30 days, with an extended tail. For the whole sample, the rise and decay timescales show a strong positive correlation, with a Spearman correlation coefficient of ρ = 0.73 and a null probability p < 10⁻⁸; i.e., slowly rising events tend to decay slowly. Applying a linear fit, we find t_decay,1/e = (1.47 ± 0.07) t_rise,1/e + (0.35 ± 2.50) days, which is similar to what was found in previous studies (Nicholl et al. 2015; De Cia et al.2018).

Zoom In Zoom Out Reset image size

It is interesting to note that four of the five SLSNe-Ib in the Yan et al. (2020) sample—namely SN 2018kyt, SN 2019hge, ZTF19acgjpgh (SN 2019unb), and ZTF20ablkuio (SN 2020qef)—have rise timescales that are 9–32 days longer than that of the rise–decay time correlation. The last SLSN-Ib, ZTF19aamhhiz (SN 2019kws), in the Yan et al. (2020) sample has a short rise time, but a much longer decay time. As noted in Yan et al. (2020), these SLSNe-Ib may have He-rich CSM, and CSM interaction may affect the LC evolution before or after the peak. Detailed modeling and discussions of the LCs are given in Paper II.

Although SNe Ic and SLSNe-I have similar postpeak spectra, they have very different peak luminosities and timescales. Figure 6 compares the peak M_g and rise times between normal SNe Ic, broad-lined SNe Ic (SNe Ic-BL), and SLSNe-I. The SLSNe-I are from our sample, while the normal SN Ic sample and the SN Ic-BL sample are taken from Barbarino et al. (2021) and Taddia et al. (2019), respectively. We exclude two bright SNe Ic from the literature, iPTF12gty and iPTF15eov, which can be better classified as SLSNe-I, as suggested by their authors. Both SN Ic samples are obtained from the (intermediate) PTF, which has a similar observation depth to ZTF. The overall distribution in Figure 6 is similar to that shown in De Cia et al. (2018). It should be noted that corrections for observational selection biases and volumetric corrections are required to interpret the distribution correctly. The peak M_g of the SLSNe-I is about 4 and 3 magnitudes brighter than those of the normal SNe Ic and SNe Ic-BL, respectively. All SLSNe-Ib have low luminosities and moderate rise times compared with normal SLSNe-I. SLSNe-I have significantly longer rise times and wider dispersions compared to those of SNe Ic. First, this is physically related to the fact that they have larger ejecta masses and more massive progenitor stars (see Paper II). Second, SNe Ic are primarily powered by radioactive decay, which has a constant timescale of energy injection. In the magnetar model proposed for SLSNe-I, the timescale of energy injection depends on the magnetic field B (∝B⁻²) and spin period P (∝P²). The LC width expands with the decrease of B (see Figure 2 of Kasen & Bildsten 2010). The large spread of their rise times could suggest that the magnetic field strengths of the neutron stars can vary over a wide range as well.

Zoom In Zoom Out Reset image size

Understanding the variation and uniformity of the SLSN-I SEDs as a function of time should shed light on the physical nature of this population of stellar explosions. Although the complete characterization of the SLSN-I SEDs is beyond the scope of this paper, we can infer some basic properties by examining the distributions of the (g − r) colors and blackbody temperatures, since both parameters are determined by the transient SEDs.

Figure 7 shows the observed color evolution as a function of time. All the colors have been corrected for reddening, but not for K-corrections, since that requires better spectral coverage than we have. The overall (g − r) color trend evolves from blue (i.e., ∼−0.3 mag, hotter temperatures) at early phases to red, and reaches (g − r) ∼ + 0.8 mag at about two months after the peak. In our sample, four events have a peculiar color (temperature) evolution, and do not follow the general trend. We discuss these outliers below. We further examine the distribution of the observed (g − r) color obtained at the maximum light in Figure 8. At the peak phases, the median observed color is close to zero, ${(g-r)}_{\mathrm{med}}=-{0.03}_{-0.11}^{+0.12}$ mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 7 and 8 show large scatters. The observed colors are not easy to interpret, because they sample different parts of the SEDs, depending on the transient redshifts. The proper measurements are the rest-frame color tracks, but reliable results would require many more spectra than our sample has. Here, we compute only the rest-frame (g − r) colors at the peak phase, using the color corrections computed from the near-peak spectra. The rest-frame (g − r) colors are corrected from the observed (g − r) for the events at z ≤ 0.17 and (i − r) for z > 0.17, if photometry data and spectra are available. All the peak colors are listed in Table A5. The green histogram in Figure 8 shows that the median rest-frame (g − r)_rest,med is $-{0.21}_{-0.12}^{+0.19}$ mag, which is consistent with the −0.27 mag calculated from the PTF SLSN-I sample (De Cia et al. 2018). The large scatter in the peak rest-frame (g − r) indicates that the SEDs of SLSNe-I may show a diverse shape, and hence a wide range of blackbody temperatures.

Figure 9 illustrates a moderate correlation (ρ = 0.52, p < 10⁻³) between the peak rest-frame (g − r) colors and the g-band absolute peak magnitudes M_g,peak; i.e., brighter SLSNe-I tend to have bluer color, which was also previously found by Inserra & Smartt (2014) and De Cia et al. (2018) in smaller samples. If combined with the PTF sample from De Cia et al. (2018), the correlation becomes stronger (ρ = 0.56, p ≈ 2 × 10⁻⁵), and can be described by a linear function, M_g,peak = (11.5 ± 2.6) × (g − r) − (18.9 ± 0.5) mag, with a 1σ error of 1.7 mag. This correlation can be used in cosmological searches of SLSNe; e.g., Inserra et al. (2021). However, this is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

To fit the blackbody temperature, we adopt a modified blackbody function, defined by ${f}_{\lambda }=\max [\,0,\,1-A\,\times (1.0-\lambda /3000.0)]\times {B}_{\lambda }$, with B_λ being the Planck function and A being the scaling factor. The modified blackbody function aims to quantitatively capture the variations in the UV spectral suppression for different events, as shown by the Hubble Space Telescope UV spectra of SLSNe-I (Yan et al. 2017, 2018). The scaling factor A is derived from the fitting at λ ≤ 3000 Å , set to zero at λ ≥ 3000 Å and fit in a range of 0–3. Larger scaling factors represent stronger suppression in UV, and A = 1 represents the SED function used in Nicholl et al. (2017). The error associated with the temperature is estimated using the Markov Chain Monte Carlo (MCMC) method.

Of the 78 events in our sample, only 15 have at least three epochs of Swift UV photometry for properly computing the blackbody temperatures. In Figure 10, the left panel shows the temperature evolution tracks for 12 events, and the right panel shows the other three with peculiar color/temperature evolutions. Although the statistics is not large, Figure 10 indicates the large temperature spread at any given phase, especially at the prepeak and peak phases. In the left panel, some SLSNe-I have high temperatures at ∼15,000 K at t ∼ −10 days, and cool down to ∼9000 K at +20 days. Interestingly, there are also two cooler SLSNe-I (i.e., SN 2019hge and SN 2019unb), with peak temperatures less than ∼10,000 K. Consistently, neither of them show clear O ii absorption lines in their classification spectra, which require a high ionization temperature (i.e., T ∼ 15,000 K; Quimby et al.2018). We perform a polynomial fit to the temperature tracks and derive T = 0.090201 t³ − 3.4306 t² − 180.30 t + 13395 K and T = − 0.052720 t³ + 0.3319 t² − 15.84 t + 9007 K for the high- and low-temperature tracks, respectively. These are shown as the blue and red shaded regions with the ±1σ uncertainty in Figure 10. Consistent with the temperature evolution, the (g − r) colors of low-temperature events are ∼0.2–0.3 mag redder than those of high-temperature ones at the peak, but become indistinguishable from those of the full sample +15 days after the peak.

Zoom In Zoom Out Reset image size

Moreover, the color distribution plots suggest that there are more low-temperature SLSNe-I. For instance, 10 of 35 (29%) events are found to have redder peak rest-frame colors than the two low-temperature events (∼0.10 mag). However, these 10 events have no UV photometry available, thus no temperature measurements. It is possible that these 10 red events could also have low peak temperatures.

Combining the temperature measurements from the PS1 sample, we do see more low-temperature events. To better quantify the temperature distribution of SLSNe-I, we choose three typical epochs (i.e., t ∼ −10 days, 0 days, and +20 days relative to the rest-frame g-band peak) when sufficient data are available from the ZTF and PS1 sample. We show the temperature distribution at these epochs in Figure 11. At t ∼ −10 days, the temperature of SLSNe-I ranges from 7000 to 23,000 K for different events, while this range is 6000–20,000 K around the peak and 6000–12,000 K (excluding one special event, SN 2019szu) at t ∼ +20 days. As SLSNe-I evolve, both the median value and the spread of the temperature become smaller. This conclusion will still hold true if we expand the temporal range from −20 to +40 days, according to the temperature trend shown in Figure 10. Furthermore, the temperatures at any of the epochs in Figure 11 are found to show flat and unimodal distributions, indicating that SLSNe-I may have a continuous temperature distribution.

Zoom In Zoom Out Reset image size

We also measure the blackbody radius of the photosphere. As shown in the top panel of Figure 12, most events show a linearly expanding photosphere from explosion to 10−20 days around the peak. This is due to the recession of the photosphere being negligible during this phase. We apply a linear fit to measure the photosphere velocity, V_phot. In Paper II, we measured the velocities from the Fe ii and O ii absorption lines in the spectra. The velocity implied from the species V_ion is expected to be higher than that from V_phot, since only the line features formed at higher velocity and lying outside the photosphere can be observed. As shown in the bottom panel of Figure 12, the distribution of ΔV = V_phot − V_ion proves that most events have a negative ΔV and their V_phot is on average lower than V_ion by about 2000–3000 km s⁻¹. Two outliers, SN 2018bgv and SN 2018kyt, are found to have significantly higher V_phot than V_ion. SN 2018bgv is the fastest-evolving event in our sample, while SN 2018kyt also exhibits relatively fast-evolving behaviors. Thus, the inconsistency between V_phot and V_ion may be due to V_phot being the average velocity measured over a period of 10−20 days, while V_ion is fast-evolving and measured at a single epoch. All of the measured temperature and radius values are listed in Table A6.

Zoom In Zoom Out Reset image size

In summary, our measurements show that the SLSN-I have a wide range of temperatures, especially at early phases. While most SLSNe-I from ZTF have temperatures over 11,000 K at the peak, and cool down rapidly, there are also many lower-temperature and slowly evolving events. This indicates that the SEDs of SLSNe-I may have diverse shapes and different evolutionary tracks.

Finally, we discuss four peculiar events that have extraordinary color and temperature evolution, as shown in Figures 7 and 10.

SN 2018hti is well sampled in UV and shows notably higher temperatures and much bluer (g − r) colors compared with other events. Its LC can be well reproduced by a magnetar model (Lin et al. 2020). This event may represent some of the population with higher temperatures and bluer SEDs. However, we caution that this event has a very high galactic extinction E(B − V) = 0.4 mag. The dust extinction corrections in the UV bands are highly uncertain.

SN 2020fvm has two almost equally bright LC peaks (at phase ∼−60 and 0 days, respectively) and the longest 1/e maxima rise timescale, of 91 days. Both its color and temperature evolutions are peculiar. The (g − r) color around the two peaks shows a similar evolutionary trend, i.e., initially evolving from red to blue before the peak, and then turning red after the peak. During the phase from −75 days to +35 days, its temperature remains almost constant, at ∼10,000 K, which could be due to the absence of sampling. Such a long timescale and double-peak evolution are unexpected for a simple radiative cooling or magnetar model. Dessart (2019) examined the color evolutions for various configurations of magnetar models, and found that almost all have fairly blue colors in the early phase, with none showing a color change of ∼0.7 mag from red to blue. The ejecta interactions with an extended H-poor CSM could be a possible explanation.

ZTF19aanesgt (SN 2019cdt) has much redder color at +20 days after the peak, and it evolves much faster in comparison to the other objects in our sample. No temperature evolution is computed, due to the lack of UV data. This event is somewhat similar to SN 2018bgv, with both belonging to the fast-evolving and CSM model–favored events (see LC modeling in Paper II), though it has redder (g − r) color. Such a red color and fast LC evolution could also be consistent with a magnetar model with high kinetic energy, as shown in Figure 13 of Dessart (2019). Nevertheless, the Fe ii velocity and the kinetic energy of the ejecta derived for SN 2019cdt are ∼14,400 km s⁻¹ and 7.6 × 10⁵¹ erg, respectively, which are both higher than the average value (see Paper II).

Finally, unlike any other SLSN-I, ZTF19acfwynw (SN 2019szu) shows unusually high temperatures that rise as the luminosity declines. In the early phases, its colors are blue, consistent with having high temperatures. Although the blackbody temperatures at late time have large errors, due to contamination from increasingly strong nebular emission lines, its spectral sequence has revealed a clear excess of UV continuum emission well past the peak phase. One possible explanation is CSM interaction, which could offer an additional heating source, boosting the emission at shorter wavelength. Another possibility is that the ejecta could be ionized by the ionizing flux from a long-lived central source (e.g., a magnetar; Margutti et al. 2017). Ionization can increase the optical opacity dominated by the electron scattering and decrease the UV opacity dominated by line transitions of metals, leading to a shift in the peak of the SED from optical to UV frequencies.

It is worth noting that SN 2020fvm and SN 2019cdt are poorly fit by simple magnetar models. It is possible that CSM interaction plays significant roles in these two systems. Additional modeling and analysis are presented in Paper II.

The bolometric LC is an important indicator of the total radiative energy of a transient and sets constraints on the possible explosion mechanisms. Due to the high photospheric temperatures of SLSNe-I, UV radiation contributes to a large fraction of their emission (Yan et al. 2017, 2018). Since only a small portion of the sample has UV data, we first need to derive an empirical bolometric correction (BC) relation, which can tie g- and r-band photometry to the bolometric luminosity, with $\mathrm{log}{L}_{\mathrm{bol}}=\mathrm{log}{L}_{{gr}}+\mathrm{BC}$. ³² Here, L_gr is the sum of the g- and r-band luminosities.

Lyman et al. (2013) derived BCs for a sample of core-collapse SNe at low redshifts, using their well-observed SEDs. We adopt a similar method and apply it to our sample. The basic concept is as follows. BC is strongly influenced by the redshift and temperature (the slope of the spectra) at each epoch. As the observed (g − r) color is also largely determined by the same two parameters, a correlation between the BC and the (g − r) color is expected to be seen. This correlation can serve as the basis for constructing bolometric LCs with only optical data.

We first compute the bolometric luminosity at the epochs when both the UV and optical data are available, using the modified blackbody SED fit described in Section 5.2. The UV scaling factor A is derived by fitting the UV to optical SEDs. For the epochs without UV data, we set A = 1.55, which is the median value of the epochs with UV data. The bolometric luminosities are shown in Figure 14 (there are only 15 events with at least three-epoch UV photometry) and listed in Table A6.

Figure 13 shows the observed (g − r) color versus the BC, i.e., logL_bol/L_gr , the luminosity ratio. We apply a third-order polynomial fit by optimizing a Gaussian likelihood function and derive an empirical relation as $\mathrm{log}({L}_{\mathrm{bolo}}/{L}_{{gr}})\,=-1.093\,{x}^{3}+1.244\,{x}^{2}-0.261\,x+0.410$, where x = (g − r). Each point is weighted by the errors of both the color and the luminosity ratio. This equation can be used to compute the BC for low-z SLSNe-I with a similar (g − r) color evolution as our sample. Note that our fit is applicable for − 0.4 < (g − r) < +0.9 mag, the phase range of −74 days < phase < + 173 days and 0.06 < z < 0.57.

Zoom In Zoom Out Reset image size

As inferred from the residuals, the systematic error for the derived bolometric luminosity is about 19%. The BC uncertainty (shown as the shaded area in Figure 13) combines the errors from the MCMC estimates and the systematic error. This error is the dominant one compared to other error sources, like host galaxy reddening and redshifts. In Figure 14, the bolometric luminosities constructed with both UV and optical photometry are consistent with the ones derived from g- and r-band photometry. This illustrates the reliability of our method.

Zoom In Zoom Out Reset image size

The bolometric LCs derived from g- and r-band photometry are shown in Figure 15. The errors of both photometry and BC are combined together in quadrature to represent the error of the bolometric LCs. ZTF19aauiref (SN 2019fiy) and SN 2019aamu are excluded from this figure, due to the bad photometry quality. Figure 15 illustrates a diversity in both the peak luminosities and LC widths. Bumpy LCs are commonly seen and we perform detailed analysis in Paper II.

Zoom In Zoom Out Reset image size

The peak bolometric luminosities are tabulated in Table A4 and the distribution is shown in Figure 16. The peak bolometric luminosities show a similar distribution as the absolute peak magnitudes. These are the raw distribution functions, without any corrections for selection biases. Further studies are required to confirm such bimodal distributions. The peak bolometric luminosity spans from 3.22 × 10⁴³ to 7.80 × 10⁴⁴ erg s⁻¹, with a median value and 1σ dispersion of ${L}_{\mathrm{bolo},\mathrm{med}}\,={2.00}_{-1.44}^{+1.97}\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Our measurements are consistent with the results from the PTF and PS1 samples. Note that the DES sample tends to have lower luminosities, which could be due to their being pseudo-bolometric luminosities constructed from the trapezoidal integration of photometry. Compared to the median peak luminosity of SNe Ib/c (∼2 × 10⁴² erg s⁻¹; Prentice et al. 2016), SLSNe-I are about 100 times more luminous, indicating different energy sources and/or explosion mechanisms from normal core-collapse SNe.

Zoom In Zoom Out Reset image size