SN 2015bn: A DETAILED MULTI-WAVELENGTH VIEW OF A NEARBY SUPERLUMINOUS SUPERNOVA

Follow-up imaging observations of SN 2015bn were triggered with a number of instruments. Optical data in ugriz were obtained with the 2.0 m Liverpool Telescope (LT) using the IO:O imager, and the 1.0 m Las Cumbres Observatory Global Telescope network (LCOGT). After ∼200 days, deep ugriz images were obtained with EFOSC2 on the NTT. These images were reduced by respective facility pipelines (the PESSTO pipeline is described by Smartt et al. 2015), including de-biasing and flat-fielding. PS1 data obtained in the ${w}_{{\rm{P1}}}$-band were converted to the more standard r-band by a shift of $+0.1$ magnitudes. The same shift was applied to CRTS R-band data to match the better-sampled r-band. These shifts were calculated from the observed spectra. UBVRI data were obtained from the 16'' Ritchey–Chretien telescope (SLT) at Lulin Observatory, Taiwan. These were reduced using standard techniques in iraf, correcting for bias, flat-field and dark current. Near-infrared (NIR) images were taken using SOFI on the NTT and NOTCam on the Nordic Optical Telescope (NOT). These were reduced using the PESSTO pipeline, and the external iraf package notcam (version 2.5)³⁰ , applying flat-fielding and sky-subtraction on dithered frames.

Photometric measurements were made on the reduced frames using a custom pyraf photometry pipeline, calling standard daophot routines to fit the point-spread function (PSF) and capture the SN flux. The nightly zero-points, in AB magnitudes, were calculated using a sequence of local field stars in the Sloan Digital Sky Survey (SDSS) Data Release 12 (Alam et al. 2015). For the SLT photometry, the SDSS ugriz magnitudes of the sequence stars were transformed into UBVRI using the equations from Jordi et al. (2006). These are reported in the Vega magnitude system. Errors in the SN magnitudes include both the scatter in the zero-points and the uncertainty in the PSF fit returned by daophot. ASAS-SN photometry was computed using a dedicated pipeline. Because SN 2015bn was close to the limiting magnitude for ASAS-SN ($\approx 17$ mag), we increased the signal-to-noise ratio of the data using a noise-weighted stacking of neighboring epochs.

For the latest epochs of ugrizJHK imaging at ∼250 days, the host galaxy contributes $\sim 20 \%$ of the observed flux. Thus we remove this flux as follows: for griz, we align the images with SDSS templates and apply image subtraction using hotpants ³¹ (based on the algorithms of Alard & Lupton 1998). For uJHK, where template images without the SN are either not available or are too shallow to effectively remove the host flux, we do not use 2D image subtraction, but instead subtract fluxes based on a model galaxy spectral energy distribution (SED) (see Section 8) from our measured magnitudes. This method works reasonably well since the host of SN 2015bn is faint and compact, and the SN does not show a large offset from the centroid.

We activated our approved Swift-GI program dedicated to the study of nearby SLSNe (PI: Margutti), as well as a number of other Target-of-Opportunity programs (PI: Inserra). We also include public data obtained under a separate program (PI: Brown). Imaging was taken with the Ultraviolet and Optical Telescope (UVOT), in the filters uvw2, uvm2, uvw1, u, b and v. The SN flux was extracted following standard methods (Brown et al. 2009) and using a 5'' aperture, and calibrated in the Swift photometric system (Poole et al. 2008). A shift of +1.02 mag was³² used to convert the Swift u zeropoint from Vega to AB system. While the transmission functions of the Swift and ground-based u-filters are not identical, the photometry obtained matches that from LT perfectly, hence we feel justified in combining these data into a single, well-sampled light curve. All photometric measurements, including epochs of observation, instruments, and reference star magnitudes, are given in the Appendix.

The observed light curves of SN 2015bn, from UV to NIR, are shown in Figure 2. SN 2015bn reached r-band maximum light on MJD 57102. At r = 16.7 mag, it exhibits the brightest apparent magnitude of any SLSN I yet discovered. The first detection from CSS is 88 days earlier, which translates to 79 days in the rest-frame at z = 0.1136. This makes SN 2015bn a slowly evolving SLSN with one of the longest rise times to date, similar to the estimate for iPTF13ehe (Yan et al. 2015). To fill in the gap from this first detection until our classification and follow-up campaign beginning at −28 days, we searched archival images from PS1 and ASAS-SN. We recovered detections in V- (ASAS-SN), i-, and z-bands (PS1), confirming the long rise. While the wavelength and temporal sampling is extremely dense during our follow-up campaign, the pre-classification photometry is relatively sparse. This makes it difficult to constrain the shape of the light curve at early times. We will address this issue in the next section.

Zoom In Zoom Out Reset image size

In Figure 3, we show the historical magnitudes measured at the location of SN 2015bn by CRTS/CSS (Drake et al. 2009). For clarity, we plot only the average detection limit, although there are many historic non-detections (the SN 2015bn host, at $R\sim 22$ mag, is well below the survey limit—see Section 8). The CSS photometry shows two possible source detections (marked as unfilled stars on Figure 3) prior to the unambiguous SN photometry (filled stars and circles). Without access to the CSS data, we cannot check whether these points are related to SN 2015bn. If real, this would be the first detection of historic variability at the location of a SLSN I (either a pre-explosion outburst or another SN in the same galaxy). We will not speculate further on the nature these early "detections," and exclude them from the rest of our analysis.

Zoom In Zoom Out Reset image size

To convert our photometry to absolute magnitudes in the rest-frame of the SN, we correct for distance modulus, extinction, and differences in rest-frame filter wavelengths (the K-correction).³³ We use the dust maps of Schlafly & Finkbeiner (2011) to correct for Milky Way extinction with $E(B-V)=0.022$, and assume negligible extinction in the SN host. This assumption is necessary, as the host emission lines are too weak to use reliably for an estimate of internal reddening (see Section 8), but also reasonable since we do not see any strong Na i D absorption from the host, which is generally thought to correlate strongly with dust extinction (Poznanski et al. 2012). K-corrections were calculated in the optical and NIR by comparing synthetic photometry on our spectra (Section 4) in the observer and rest frames. This was facilitated by the K-correction code snake (Inserra et al. 2016). At this redshift, the size of the correction was typically $\lesssim 0.2$ mag. In the NUV, where we have no spectroscopic coverage, the correction was estimated as ${K}_{{\rm{UV}}}=-2.5\;$ ${\mathrm{log}}_{10}(1+z)\;=-0.12$ mag.

SN 2015bn has ${M}_{g}=-22.0\pm 0.08$ mag (AB) and ${M}_{U}=-23.07\pm 0.09$ mag (Vega) at maximum light (where the uncertainty includes a systematic error dominated by the K-correction as well as the statistical error), making it comparable to some of the most luminous SLSNe I (Vreeswijk et al. 2014; Nicholl et al. 2015b). However, this is about 1.5 mag fainter than ASASSN-15lh, the brightest known SN (Dong et al. 2016). The rest-frame optical light curves are shown in Figure 4. Overlaid for comparison is another well-observed, low-redshift SLSN I: PTF12dam (Nicholl et al. 2013). The decline rate after $\approx 50\;{\rm{days}}$ is similar between these two objects, but the light curves do have a number of differences. Most obvious is the slower rise exhibited by SN 2015bn, where the broad peak looks more symmetrical around maximum. SN 2015bn is also brighter than PTF12dam by $\approx 0.5$ mag in the red filters, but the peak luminosities match in u band. Thus SN 2015bn is brighter but redder than PTF12dam. However, after ∼100 days from peak, SN 2015bn is brighter in the g-band (PTF12dam was not observed in u-band at this epoch), demonstrating a very slow color evolution.

Zoom In Zoom Out Reset image size

The slow color evolution is shown explicitly in Figure 5, and compared to other SLSNe (all objects have been K-corrected to rest-frame using their observed spectra). The color at peak is similar for most objects, with $-0.1\gtrsim g-r\gtrsim -0.4$ mag (though SN 2011ke seems to have $g-r\gtrsim 0$ mag). However all the objects with slowly declining light curves remain bluer for much longer than their fast-declining counterparts. SN 2011ke, which settled onto a slow light curve tail after around 50 days from maximum light, has a comparable g − r color to the slowly declining objects at 200 days. The u − g color of SN 2015bn initially evolves to the blue, reaching a minimum of $\sim -0.1$ mag at maximum light, before declining relatively steeply compared to the optical-only colors. The faster decline in the UV is expected as the SN cools. Optical colors are roughly constant from our first observations at −25 days until shortly after peak, with $g-r\approx r-i\approx i-z\approx -0.2$ mag. The optical colors then slowly evolve to the red. SN 2015bn reaches $g-r\approx 0.4$ mag by +250 days. This evolution is slower and bluer than any literature SLSN with such late observations, as is seen in the bottom panel of Figure 5, and from the g-band excess compared to PTF12dam in Figure 4.

Zoom In Zoom Out Reset image size

Many SLSNe I exhibit a non-monotonic rise to maximum brightness, with a fainter initial peak, or "bump." These bumps have been clearly detected for a number of individual objects (Leloudas et al. 2012; Nicholl et al. 2015a; Smith et al. 2016). Recently, Nicholl & Smartt (2016) have shown that a number of literature objects are also consistent with an undersampled bump, and the relative faintness of this first peak makes it difficult to exclude in most cases. In Figure 6, we plot the rising part of the light curve in the two filters with the best early sampling: r- and V-bands. We fitted a third-order polynomial (top panel) to each light curve. While either band alone would be insufficient to test for the presence of a bump, the two bands together show that if we were to assume a smooth rise, we would derive a very surprising color evolution. The V-band data clearly show a steeper decline than that implied by a smooth fit to r-band. The V − r color between −30 days and peak is constant at −0.15 mag, but the fits suggest the color at early times would be more like $V-r\approx 1$ mag. Such a dramatic change in color over the rising phase (and such a red color at early times) is highly inconsistent with other SLSNe I that have multicolor data available at this phase (Leloudas et al. 2012; Howell et al. 2013; Inserra et al. 2013; Smith et al. 2016). Such objects generally show blue colors, similar to what we see around the light curve peak here, and a flat color evolution before maximum. By contrast, if we assume that the color evolution is relatively flat, as in other SLSNe I, and combine our V- and r-band data into a single light curve using the observed color at peak, we find that the rising light curve of LSQ14bdq (Nicholl et al. 2015a) gives an excellent match to the data, after applying a temporal stretch factor of 1.2 (bottom panel). Thus the multicolor data seem to be more consistent with a bump similar to LSQ14bdq and other SLSNe I observed at very early times. However, we note that this is a larger stretch than was required for any of the objects in Nicholl & Smartt (2016).

Zoom In Zoom Out Reset image size

Our high-cadence photometry reveals the presence of several more distinctive bumps/plateaus in the light curve, most prominent in the UV bands. First the u-band shows a plateau lasting 14 days in the rest-frame, until 10 days before maximum light (Figures 2 and 4). This shows up clearly as a dip in the UVOT uvm2 light curve, which has an effective wavelength ∼2230 Å in the observer frame, or ∼2000 Å at the redshift of SN 2015bn. We designate this feature "the shoulder" for all subsequent discussion. A second undulation is present in all optical and UV filters at 50 days after maximum light, or MJD 57158—we call this "the knee." It is known that UVOT observations can occasionally be affected by artificial drops in flux (Edelson et al. 2015). However the variability exhibited by SN 2015bn is reliable for a number of reasons. These flux "drop-outs" are not expected to show correlated deviations spanning numerous successive epochs. Furthermore, the biggest deviation is seen in uvm2, which was found by Edelson et al. (2015) to have a lower drop-out rate than uvw2 or uvw1. Most importantly, these light curve fluctuations are seen at the same phases not only across all of the UVOT filters, but also match perfectly (in time and brightness) with the fluctuations in our ground-based data.

To see the undulations more clearly, we fitted third-order polynomials to the rest-frame light curve decline in each UV and optical filter (after masking epochs between +40 and +60 days), and subtracted the fit from the full declining light curve. Unfortunately, the NIR cadence is too sparse to show if the knee exists there too. The residuals are shown in Figure 7. The excess is more pronounced in the blue and UV bands. This temporary evolution to the blue can also be seen in Figure 5, where the u − g color evolves from ∼0.6 mag at +30 days, to ∼0.4 mag at +50 days. A third undulation may occur at +100 days. This is quite prominent in g-band, but the UV light curves are not well sampled at this epoch. This type of behavior—a temporary late-time re-brightening, more pronounced in the blue—has been seen before in one other SLSN, SSS120810 (Nicholl et al. 2014). SN 2007bi also seemed to exhibit a bump at $\approx 100\;{\rm{days}}$ after maximum, but no color information is available at this phase (Gal-Yam et al. 2009). The implications of these fluctuations will be discussed in Sections 5 and 6.

Zoom In Zoom Out Reset image size

Optical spectra of SN 2015bn were acquired using the NTT and EFOSC2, the 6.5 m Magellan telescopes with the IMACS and LDSS3 spectrographs, WiFeS on the Australian National University 2.3 m telescope (Dopita et al. 2007), LT with SPRAT, the 2.0 m LCOGT Faulkes Telescopes and FLOYDS, and the University of Hawaii 2.2 m telescope with SNIFS. All spectra were reduced, including de-biasing, flat-fielding, object extraction, wavelength and flux calibration, using instrument-specific pipelines or standard routines in iraf. The relative flux calibrations were achieved using standard star observations taken on the same nights as the SN spectra. Details of the complete list of spectroscopic observations are given in the Appendix. NIR spectra were obtained using SOFI on the NTT, and reduced using the PESSTO pipeline to flat-field, sky-subtract and co-add the frames, apply wavelength calibration to the 2D frames, as well as a flux calibration and telluric correction from standard star observations, and finally to extract the 1D spectrum. For every spectrum, the absolute flux calibration was checked against contemporaneous multicolor photometry; spectra were then scaled to match the photometry if necessary. PESSTO raw data are immediately available in the ESO archive³⁴ , and the reduced data for the season in which SN 2015bn was observed are due for bulk public release in late 2016.³⁵ In the meantime, all spectra are available on WISeREP.

The complete set of optical spectra are shown chronologically in Figure 8. SN 2015bn retains a blue continuum throughout the observing season, and this is still present as late as 240 days after maximum light, i.e., the ejecta are not yet fully nebular. The lines in the spectrum also evolve very slowly. To show more clearly how the individual lines evolve, and to identify the ions responsible, we normalized each spectrum by dividing by its mean flux, and fitted the continuum with a third-order polynomial, which was then subtracted off. Plotting the normalized, continuum-subtracted spectra in Figure 9, we can see directly how the line centers, velocities and strengths evolve with time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The line identifications in this section were made through detailed comparison with identifications and SYNOW fits (Thomas et al. 2011) given for other SLSNe I by Gal-Yam et al. (2009), Young et al. (2010), Pastorello et al. (2010), Quimby et al. (2011), Chomiuk et al. (2011), Leloudas et al. (2012), Inserra et al. (2013), Lunnan et al. (2013), Nicholl et al. (2013, 2014). Furthermore, we examined similar modeling of normal SN Ic spectra by Millard et al. (1999), Sauer et al. (2006), Valenti et al. (2008), Hunter et al. (2009). Line wavelengths in the rest frame were determined by consulting the NIST Atomic Spectra Database (Kramida et al. 2015), and compared to absorption lines assuming a photospheric velocity of 8000 km s⁻¹ (see Section 4.3). For the most common ions, multiple lines of the same species were checked for consistency using the interactive plotting functionality of WISeREP (Yaron & Gal-Yam 2012). A full set of self-consistent models to reproduce the spectral evolution of SN 2015bn will be the subject of future work, but is beyond the scope of this paper. In the discussion below, we use the symbol λ when referring to lines by their rest-frame wavelength, and Å when referring to redshifted absorption.

It is immediately clear from Figures 8 and 9 that the spectra divide fairly cleanly into two phases. Very little evolution is seen from our first spectrum at −27 days through to +7 days. Our highest signal-to-noise spectra taken over this period are overlaid in Figure 9 (purple to dark green lines). Between +7 and +20 days, the spectrum undergoes something of a transformation, after which it evolves very slowly for the next 200 days. This result is independent of any model assumptions or line identifications. Sample spectra from this second phase are also overlaid in the figure (light green to red lines).

The early ($t\leqslant 7$ day) spectra are dominated by blue continuum emission, superimposed with fairly weak, broad absorption features and P Cygni lines. In the blue, we clearly see the O ii absorption lines that have come to characterize the SLSN I class; however, in this case the absorption at 4300 Å is stronger and broader than the other lines. This could indicate a contribution from Fe iii λ4430. Between 3000 and 3600 Å, we detect absorptions that have been attributed to Fe iii and Si iii in other SLSNe I (e.g., Leloudas et al. 2012; Lunnan et al. 2013). Moving to redder wavelengths, we see signatures of Fe ii at ∼5000–6000 Å, along with Na i D and Si ii λ6355. We tentatively associate the feature at ∼7000 Å with C ii λ7234. To the red of this we see a fairly strong O i λ7774, which may include some contribution from a Mg ii line of similar wavelength.

We may also see some spectroscopic indication of hydrogen at early times. Yan et al. (2015) presented an Hα emission line in the nebular spectrum of iPTF13ehe, caused by late-time ($\gtrsim 200$ day) interaction with hydrogen-rich material at a large radius, but no hydrogen was seen in the photospheric phase. The upper right panel of Figure 9 includes the region around Hα (rest-wavelength 6563 Å). There is an absorption in the earliest spectrum at 6300 Å, which is consistent with Hα at a velocity $\approx 8000$ km s⁻¹. However, an alternative identification for this line is C ii λ6580, particularly if the line at 7000 Å during this phase is also carbon. On the other hand, the absorption at 6300 Å vanishes by the time of maximum, light, whereas the 7000 Å line persists.

In a further attempt to break the degeneracy, we look more closely at the region around Hα and Hβ (in velocity space) in Figure 10. The strongest line in the vicinity of Hβ is Fe ii λ4924, which was identified for SN 2007bi by Gal-Yam et al. (2009), Young et al. (2010), and is consistent with a SYNOW fit to LSQ12dlf by Nicholl et al. (2014). The initial blueshift of the P Cygni peak may indicate that the line is at first blended with Hβ; this blueshift quickly disappears along with Hα. There is also an absorption in the earliest spectrum located at $\approx 8000$ km s⁻¹ bluewards of Hβ. However, the effects of line blending make a firm association with Hβ difficult. Overall, we consider it plausible, though by no means certain, that these lines are from unburned hydrogen. Parrent et al. (2015) showed that the line profiles around 6200-6500 Å in superluminous and normal SNe Ic may be better reproduced by including hydrogen as well as Si ii. Although their analysis was based on a SLSN spectrum well after maximum light, the resolved absorption features in the early spectrum of SN 2015bn may provide independent evidence of this claim.

Zoom In Zoom Out Reset image size

If hydrogen is present in the spectrum of SN 2015bn, indicating that some small amount of the progenitor's envelope managed to remain bound before explosion, we might also expect to see some helium (though this is much more difficult to excite). While the strongest optical He i line is generally hard to distinguish from the Na i D lines, our NIR spectra (Figure 11) show a clear line at $\approx 10800$ Å that could be He i λ10830. This line is expected to be 2–10 times stronger than the strongest optical line of He i (Inserra et al. 2013). Few NIR spectra of SLSNe I exist, but Inserra et al. (2013) saw a similar feature in SN 2012il. Normal SNe Ic generally exhibit P Cygni lines with a deep absorption component at this wavelength. However, the identification of this line has been hotly debated in the literature. (Patat et al. 2001) claimed a detection of He in the NIR spectrum of the LGRB-associated Type Ic SN 1998bw, but this has been challenged by Taubenberger et al. (2006), Hachinger (2011). It has been argued that this line may contain a large contribution from C i, since the other strong helium line in the NIR ($\sim 2.1\;\mu {\rm{m}}$) tends to be weak or absent in SNe Ic, unlike in He-rich SNe Ib (Hunter et al. 2009). Unfortunately, our wavelength coverage does not extend to 2 μm for SN 2015bn. While the line could plausibly be He i, a firm identification is not possible. The presence of carbon in the optical spectrum may suggest that C i is a more likely explanation for this feature. Other species that could contribute at this wavelength are Mg ii and Si ii. Future spectral modeling will be required to determine the nature of the 10800 Å line.

Zoom In Zoom Out Reset image size

Returning from our digression to the NIR, we now examine the optical evolution after maximum light. Looking again at Figure 9, the most obvious changes occurring between 7 and 20 days after maximum are the emergence of Mg i] λ4571 emission, strong Fe ii lines between 4000 and 5500 Å, and Na i D absorption. There is also an unidentified line at around 6100 Å—this could perhaps be a detached, high-velocity component of Si ii with a velocity of ∼16000 km s⁻¹, but it would be very surprising for a high-velocity line to appear only at late times. At the same time, the O ii and Fe iii lines that previously dominated the blue part of the spectrum weaken and disappear as the ejecta cool and these species recombine. The P Cygni line at 7000 Å is replaced by an asymmetrical emission line due to [Ca ii]λλ7291, 7323. The [Ca ii] shows a steep blue side but a broad red wing. One explanation for this asymmetry could be blending with C ii or Fe ii. At the same time, we start to see likely emission from [O i]λλ6300, 6364 and the Ca ii NIR triplet at 8500 Å. The strengths of the calcium lines, as well as O i λ7774 and the Na i D absorption, increase substantially after +50 days. The final spectrum obtained at +243 days (and after SN 2015bn had returned from behind the Sun) remains largely similar, with significant continuum persisting in the blue. The only noticeable change over this period is that some emission lines (especially calcium) increase in luminosity, and absorption lines (especially sodium) increase in depth. The spectra at $\gt 100$–200 days, with necessarily longer exposure times, finally reveal the presence of faint nebular emission lines originating in the host galaxy, from which we securely derive the redshift of z = 0.1136.

The various subplots in Figure 9 show close-ups of different spectral lines. We can see that the minima of all absorption lines, as well as the peaks of some P Cygni profiles, move gradually to lower velocities (i.e., redder wavelength) as the spectrum evolves. We measured the photospheric velocity as a function of time using a number of ions by fitting Gaussian profiles to each line and determining the blueshift of the absorption minima. It is common to use the Fe ii λ5169 line to measure photospheric velocities (e.g., Nicholl et al. 2015b). However, in SN 2015bn this line seems to be blended with Fe iii at early times. Moreover, our densely sampled spectra, spanning ∼200 days in the rest-frame, shows that the peak of the P Cygni profile moves significantly to the red as the spectrum evolves, such that we cannot be confident that the same line dominates the absorption at all times. To get the most consistent velocity measurements, we choose lines that appear to be unblended for the majority of spectra. The smooth evolution toward lower velocity is particularly clear for the Si ii λ6355 and O i λ7774 lines, so we used these two lines, along with Fe ii λ4924, to measure the photospheric velocity of SN 2015bn.

The velocity measurements are shown in Figure 12, and compared with other SLSNe (using the same lines). We performed the measurements ourselves for PTF12dam and SN 2010gx, whereas the velocities for SNe 2007bi and 2011ke were taken from Young et al. (2010) and Inserra et al. (2013) respectively. All of the ions indicate relatively lower velocities in SN 2015bn: $\approx 8000$ km s⁻¹ at maximum light compared to 9000–12,000 km s⁻¹ in the other objects. Nicholl et al. (2015b) estimated a median maximum-light velocity of 10,500 ±3000 km s⁻¹ for SLSNe I, so SN 2015bn falls at the low-end of the expected range. The O i line shows the highest velocity, and Fe ii the lowest, as we would expect for a sensible stratification of the ejecta, with heavier nuclei toward the center. The effect is compounded as the O i line is relatively strong, such that absorption in the outermost, rarified material (which has the highest velocity) is still significant.

Zoom In Zoom Out Reset image size

The temporal evolution in the measured velocities is very slow, decreasing by only ∼3000 km s⁻¹ during a period of nearly 150 days, and seemingly staying constant thereafter. For comparison, the Fe ii velocity of the normal Type Ic SN 2007gr is also shown. Although these SNe show similar post-maximum velocity, SN 2007gr underwent a steep deceleration around maximum light. The relatively flat velocity curves of SLSNe I have been interpreted as one possible signature of a central engine accelerating the inner ejecta (Kasen & Bildsten 2010; Chomiuk et al. 2011; Nicholl et al. 2015b). The unprecedented spectroscopic coverage of SN 2015bn—together with low velocities that reduce some of the effects of line blending—reveals the clearest indication yet that at least some SLSNe have very shallow velocity gradients in the ejecta.

In Figures 13 and 14, we compare our spectra of SN 2015bn to other SLSNe from the literature. First we consider objects with broad light curves similar to SN 2015bn (Figure 4). SN 2007bi (Gal-Yam et al. 2009; Young et al. 2010) and PTF12dam (Nicholl et al. 2013) are some of the best observed examples. As demonstrated by Figure 13, the spectroscopic evolution of SN 2015bn is an excellent match to these objects. In particular, all three objects exhibit a distinctive "trio" of lines between ∼5700–6500 Å at late times. However, there are some differences. SN 2007bi shows stronger emission lines at ∼100 days, especially in Mg i] and [Ca ii], suggesting that SN 2007bi develops a significant nebular component earlier than the other two objects.

Zoom In Zoom Out Reset image size

SN 2015bn also exhibits a slight difference compared to PTF12dam at very early times. While the two objects look virtually identical at peak, the shapes and equivalent widths of the O ii lines at 4000 and 4600 Å in PTF12dam at −15 days were very similar to those in fast-declining SLSNe I, whereas in SN 2015bn these lines are weaker and may be blended with Fe iii λ4430 at the time of our earliest spectra $\gt 20\;{\rm{days}}$ before maximum. Given that the continuum temperature looks virtually identical to PTF12dam, the differences may be due to velocity structure or composition. That the spectrum of SN 2015bn looks unchanged for longer around maximum light than PTF12dam is consistent with its shallower rise, and may be linked to the pre-peak shoulder in the light curve. In summary, the early detection and high-cadence spectra reveal a small degree of diversity among these objects, but overall the spectroscopic evolution of SN 2015bn confirms the similarity of this object to PTF12dam and SN 2007bi. Given the relatively low velocities and high signal-to-noise spectra of SN 2015bn, this is an ideal data set for detailed spectroscopic modeling in the future, to confirm our line identifications and explain the slight differences between the SLSNe. However, such modeling is beyond the scope of this paper.

We now compare and contrast the same sequence of spectra with some typical fast-declining SLSNe I (the same objects used for the color comparison in Figure 5). In fact, the most obvious difference in Figure 14 was already revealed by Figure 5: that while most SLSNe I evolve to the red after maximum light, SN 2015bn (and other slow decliners) maintain a strong blue continuum for much longer. The other important difference is in the line widths; the fast-declining objects show significantly more Doppler-broadening and blending compared to SN 2015bn, as expected from their higher velocities (Figure 12). Despite these differences, the slow- and fast- declining SLSNe I do display nearly all the same lines, and in the same order—particularly O ii before maximum light and Ca ii, Mg i], Fe ii and Si ii afterwards—but the lines only become prominent at $+30\mbox{--}50\;{\rm{days}}$ in objects like SN 2015bn, in keeping with their slower photometric evolution. This degree of similarity clearly suggests a link between the fast and slow SLSNe I. Higher ejecta mass could simultaneously explain the broader light curve, lower velocity, and later formation of lines in the slow objects like SN 2015bn. In order to account for the persistent blue continuum, more late-time heating may be required.

Zoom In Zoom Out Reset image size

The literature contains numerous comparisons between the early (photospheric) spectra of SLSNe I and of normal and broad-lined SNe Ic (e.g., Pastorello et al. 2010; Inserra et al. 2013). Given the high signal-to-noise late-time spectra of SN 2015bn, we now have an opportunity to do a similar comparison times when the objects are evolving toward the nebular phase. A full nebular comparison is not possible yet, as SN 2015bn retains significant continuum and absorption lines even at +243 days. In Figure 15, we plot the +106 day spectrum of SN 2015bn along with the spectroscopically normal (but high mass) Type Ic SN 2011bm (Valenti et al. 2012) and the broad-lined, LGRB-associated Type Ic SN 1998bw (Patat et al. 2001) at comparable phases. We also show the spectrum of the fast-declining SLSN I SSS120810 at +44 days. Unfortunately, no high signal-to-noise spectra of this kind of event is available beyond $\approx +50\;{\rm{days}}$, as high redshifts make observations of fast-declining SLSNe challenging at this phase. A later spectrum of SSS120810 at +60 days does exist (Nicholl et al. 2014) and looks virtually identical to this one, but we use the earlier spectrum as it has a much higher signal-to-noise ratio. Given the faster spectroscopic (and photometric) evolution of these events, it is not unreasonable to compare a fast event at $\approx +50\;{\rm{days}}$ with SN 2015bn at $\approx +100\;{\rm{days}}$.

Zoom In Zoom Out Reset image size

The four objects show all the same major spectral features, apart from in the region between ∼5700–6500 Å, where SN 2015bn shows the distinctive trio of lines. The fast-declining SLSN does not show Na i D or the unidentified line at this phase. We note that SN 2015bn started to exhibit the trio of lines as early as +50 days, despite the much slower evolution in all of the lines it shows in common with SSS120810. This gives us some confidence that this trio may represent a real spectroscopic difference rather than an artefact of the phases we choose for the comparison. The GRB-SN 1998bw does not exhibit the central, unidentified line of the trio, but agrees in the other two lines (Si ii and Na i). This could be a consequence of line blending due to the higher velocities in this object. SN 2011bm does show three absorptions within the shaded region, but apart from Na i, the profiles seem to be redshifted relative to SN 2015bn, and may not actually be the same lines. The other distinctive difference is in the Ca ii lines. SN 2015bn is the only object that shows a larger luminosity in the forbidden line at 7300 Å than in the allowed NIR triplet at 8600 Å. A detailed explanation of these line differences is beyond the scope of our study, but may hold clues to whether there are significant differences in the ejecta conditions between different stripped SNe. In particular, the strong [Ca ii] emission in all slowly declining SLSNe has yet to be understood (Gal-Yam et al. 2009).

The next step in our analysis was to construct the bolometric light curve of SN 2015bn. Because the photometric observations span the full wavelength range from the UV to the NIR, the bolometric luminosity could be estimated relatively straightforwardly. After applying reddening corrections and K-corrections to the measured magnitudes, and correcting for distance modulus at z = 0.1136, the absolute rest-frame light curve in each filter was interpolated to match the epochs with g- and V-band observations. Magnitudes were then converted to spectral luminosity (L_λ) to give a SED at each point on the light curve. To get the luminosity, we integrated each SED numerically, assuming fluxes go to zero at the blue edge of the uvw2 band and the red edge of K-band.

We also fitted blackbody curves to each SED. As well as enabling a check that the SED looks sensible, these fits were used to derive temperatures and radii. In Figure 16, we plot the SED and blackbody fits at 4 representative epochs: our earliest epoch with multicolor photometry at −27 days; maximum light; +30 days, sampling the steeper decline prior to the knee; and +75 days, sampling the slower post-knee decline. At early times (top row of figure), the blackbody fits are poor around the peak of the SED. This has been observed before (e.g., Chomiuk et al. 2011), and is a consequence of metal line absorption in the UV, which significantly suppresses the flux bluewards of ∼3000 Å. Therefore we also fitted a "two-temperature" blackbody model, using two components that meet at the u-band. These curves (shown in Figure 16) are not intended to suggest two physically distinct temperature zones, but simply to demonstrate the very different color temperature in the UV due to the strong absorption. The red component (u- to K-band), is not subject to substantial line blocking, and should thus give a better overall representation of the color temperature of the underlying continuum emission. The amount of line absorption increases with time across the UV and optical (Figure 9). However, as the ejecta cool and the peak of the SED moves into the optical, we find that by ∼75 days the SED can be described well by a single blackbody fit from the UV to the NIR. At all epochs, integrating the two-component blackbody model from 1000 to 25,000 Å gives a luminosity estimate that is consistent (to within $\lt 5$%) with that obtained by directly integrating the observed fluxes.

Zoom In Zoom Out Reset image size

In the UV and optical filters, our photometric sampling is sufficient that linear interpolation between epochs was generally only over periods of a few days, and therefore does not wash out the substructure in the light curve (shoulder and knee; see Section 3). The NIR is not sampled quite so densely, but this is not problematic, as the NIR contribution to the overall flux is relatively modest. As demonstrated by Figure 16, at maximum light the NIR contribution is $\lesssim 10 \%$, rising to $\sim 20 \%$ during the decline phase. Since we lack NIR observations earlier than −9 days, we assumed constant colors in $g-J,g-H$ and g − K for the early points. The blackbody fit in the first panel of Figure 16 shows that this is a reasonable assumption. At the very earliest epochs, we have only a single filter—in this case we assumed the same bolometric correction as at −27 days. This assumption should be perfectly reasonable for small extrapolations, but is questionable for the earliest CRTS point—for example, Nicholl & Smartt (2016) and Smith et al. (2016) found temperatures $\approx {\rm{25,000}}$ K during the early bump phases of SLSNe, compared to temperatures of 10,000–15,000 K at maximum light. For our latest photometric point at +243 days, we lack UV data. The bolometric luminosity here was estimated by integrating a one-temperature blackbody fit to the optical and NIR data.

The bolometric light curve, as well as the temperatures and radii inferred from the blackbody fits, are plotted in Figure 17. We divide the evolution into five phases for our analysis; these are labelled a, b, c, d and e in the figure. In region a, the evolution is only partially observed. We have an isolated point at −79 days with a large error due to the reliance on an uncertain bolometric correction. At this phase, the point could be part of a long smooth rise, but our analysis in Section 3 suggested that it was more likely to be connected with the initial bump phase exhibited shortly after explosion by many SLSNe I (Leloudas et al. 2012; Nicholl et al. 2015a; Nicholl & Smartt 2016; Smith et al. 2016). When we pick up the object again with ASAS-SN, it is −54 days from peak and rising. As we only have single-filter photometry at this time, we have no information about the temperature or radius during region a.

Zoom In Zoom Out Reset image size

Region b encompasses the pre-maximum shoulder, first pointed out in Section 3. This lasts for around 12 days, and seems to be a flat plateau in bolometric space. The blackbody temperature decreases from ∼12,000 to ∼11,000 K during this phase (based on optical-NIR blackbody fit). This drop is compensated by an increase in the radius of the emitting surface. The change in radius over time during this phase is consistent with our measured velocities in Figure 12—these curves are overplotted here. If we extrapolate the O i velocity backwards to ${R}_{{\rm{BB}}}=0$ (assuming a deceleration of 26 km s⁻¹ day⁻¹, based on the velocity evolution shown in Figure 12), we arrive at an explosion date 92 days before maximum light. However, as the ejecta will likely have undergone steeper deceleration at early times, the 92 day rise should be considered an upper limit only. If the explosion occurred $\lt 92\;{\rm{days}}$ before peak, then the detection at −78 days would most likely be during the initial bump phase, assuming SN 2015bn follows a similar morphology to LSQ14bdq. This agrees with our analysis in Section 3.

The next phase of the light curve, labelled c, is the broad, smooth peak. Over 15 days, SN 2015bn rises to a bolometric maximum with ${L}_{{\rm{peak}}}=2.3\pm 0.4\times {10}^{44}$ erg s⁻¹. The light curve is quite symmetric close to maximum, then settles onto a smooth decline until ∼30 days later. The luminosity closely tracks the temperature evolution, which shows an increase back to 12,000 K at maximum and then declines smoothly afterwards. The radius largely continues to track the line velocities, indicating that photospheric recession is relatively unimportant during this phase. Therefore the changes in the spectrum occurring a few weeks after maximum light are primarily driven by the decrease in temperature. That the spectrum shows so little evolution prior to this can also be understood in terms of the temperature evolution, which varies by only around 10% from our first measurement until ∼15 days after maximum.

The most interesting behavior occurs at $\gtrsim 30\;{\rm{days}}$. This is the beginning of the knee in the light curve, which we designate region d. The blue excess visible in Figure 7 shows up clearly as a brief plateau in the bolometric output, during which time the temperature stays constant at around 8500 K. The radius stops tracking the measured line velocities, and begins to decrease from its maximum value of 10¹⁶ cm, indicating that the photosphere is receding through the ejecta faster than the ejecta are expanding. Up to this point, the measured velocities were slowly decreasing, but during region d they settle onto a constant ∼6500 km s⁻¹ (see also Figure 12). This suggests a change in the ejecta structure, with a flatter velocity profile inside some mass coordinate. This will be discussed in more detail in Section 6.5.

Finally, in region e the bolometric luminosity settles onto a long, slow decline, which we followed for 60 days before SN 2015bn disappeared behind the Sun. Since its return just before the end of 2015, our observations are consistent with the same slow decline rate, as demonstrated by a linear fit across this part of the light curve. The temperature is approximately constant at 7000 K. Meanwhile the radius looks to have decreased slowly and smoothly to ∼4 × 10¹⁵ cm. During this phase, the blackbody fits to the full UV-NIR wavelength range give much the same color temperature as the optical-NIR fits; this is expected given Figure 16. The total radiative output integrated across regions a–e is ${E}_{{\rm{rad}}}=2.3\pm 0.5\times {10}^{51}$ erg.

In Figure 18, we compare the bolometric light curve of SN 2015bn to other SLSNe I. The late-time decline rate of SN 2015bn is a near-perfect match for PTF12dam and SN 2007bi. The fast-declining SLSN 2011ke transitioned to a tail phase after around 50 days, which also matches this slope (though not all SLSNe seem to show this behavior). The decline rate is suspiciously similar to the decay of ⁵⁶Co; this will be discussed in more detail in Section 6. Looking more closely at the light curve shape around peak, SN 2015bn does show a marked difference from other slow-declining SLSNe I observed so far. PTF12dam and SN 2007bi appear to begin their slow-declining tail phase immediately after maximum light, whereas SN 2015bn initially shows a steeper decline more reminiscent of faster SLSNe. The measured bolometric decline rate between days +10 and +30 is 0.038 mag day⁻¹, which is almost identical to LSQ12dlf around the same phase, whereas PTF12dam declines at a rate of only 0.011 mag day⁻¹. The initially steeper decline rate may indicate a lower mass compared to PTF12dam; however, the rise time may be longer, and argues in the opposite direction. Mass estimates from model fits will be presented in Section 6.

Zoom In Zoom Out Reset image size

While on the subject of the rise time, another useful comparison can be made to PTF12dam if we assume that the earliest detection of each event is close to the time of explosion; this is shown in the bottom panel of Figure 18. In this case, the time to peak of PTF12dam matches the time to the shoulder in SN 2015bn. The discrepancies between the light curves then appear between ∼60 and 140 days after explosion, as SN 2015bn undergoes a number of undulations of decreasing amplitude (there is potentially a third, weak undulation around 150 days), before joining smoothly onto a tail phase that again matches PTF12dam. This shows that, without a clear picture of the physical mechanism underlying the light curve undulations, it is difficult to know which phases are most appropriate to make like-for-like comparisons between the SNe. That the spectrum of SN 2015bn during the shoulder phase looks like PTF12dam at maximum light may indicate that this alignment gives a fairer comparison than the "days-from-peak" approach used elsewhere in this work. Fortunately, this is not crucial for the analysis to follow.

To better understand the nature of SN 2015bn, we fitted the bolometric light curve with diffusion models based on the various power sources proposed for SLSNe: a magnetar engine; circumstellar interaction; and radioactive decay. The necessary equations and fitting procedures have been extensively described by Inserra et al. (2013), Chatzopoulos et al. (2012), Nicholl et al. (2014), and are based on the original diffusion solution of Arnett (1982) with constant opacity. There is a degree of uncertainty as to what value to use for the opacity when modeling SLSNe. In the optical, the dominant source of opacity is electron scattering. For normal SNe Ic, most authors assume $\kappa =0.1$ cm² g⁻¹. In the case of SLSNe, ionization may be more complete, especially before and around maximum when the temperature is higher than in normal SNe Ic; full ionization gives $\kappa =0.2$ cm² g⁻¹. In the present work we will assume this latter value. For fixed kinetic energy, the value of ${M}_{{\rm{ej}}}$ derived from fitting the light curve scales as ${M}_{\mathrm{ej}}\propto {\kappa }^{-2/3}$. Therefore our ejecta mass estimates can easily be compared with literature models that use $\kappa =0.1$ cm² g⁻¹ simply by multiplying them by a factor of 1.6.

We first consider models in which the luminosity is powered by magnetar spin-down. This is a specific example of the more general class of "central engine" models, which also include accretion-powered SNe (Dexter & Kasen 2013). The magnetar model assumes that the core-collapse of the SLSN progenitor leaves behind a highly magnetized neutron star (B ∼ 10¹⁴ G) rotating with a period $P\sim 1\mbox{--}10$ ms. The nascent pulsar thus has a large reservoir of rotational energy, which is tapped via the B-field and heats the ejecta at the rate prescribed by the standard magnetic dipole formula (Ostriker & Gunn 1971; Kasen & Bildsten 2010; Woosley 2010). The magnetar emission is expected to be in the form of energetic electron-positron pairs, which generate a hard input spectrum through a pair-cascade (Metzger et al. 2014). We treat γ-ray leakage from the ejecta following Wang et al. (2015) (see also Chen et al. 2015).

As discussed in Sections 3 and 5, because the early light curve is not well sampled, we do not know for sure whether SN 2015bn rises monotonically from discovery to maximum light, or if the first detection is during the initial bump phase that seems to be common in SLSNe I (Nicholl & Smartt 2016). The color evolution suggested that the bump interpretation is more likely (see Section 3). The basic form of the magnetar model employed here does not accommodate a non-monotonic rise (though Kasen et al. 2015 have shown that including delayed shock breakout driven by the magnetar can give a good fit to the double-peaked light curve of LSQ14bdq). Nicholl et al. (2015b) showed that, for both magnetar models and observed SLSNe, the rise and decline timescales obey a fairly tight correlation. Thus the nature of the first detection, driving the fit at early times, will have a strong influence on the fit at peak, and therefore any derived parameters. For this reason, we calculated the best-fitting model under two different constraints: including the first point as part of a smooth rise; and introducing a synthetic point around where we would expect the monotonic rise to begin if the first detection is in the bump phase—we do this by scaling the light curve of LSQ14bdq as in Figure 2.

The model fits are shown in Figure 19, with the fit parameters listed in Table 1. In general, the fits are seen to be reasonable; however these simple models are not able to replicate the detailed structure in terms of the knee and shoulder. The value of the derived mass is driven by two competing factors. The relatively narrow shape around the main light curve peak can be more easily fit with a lower-mass model, while the overall broader shape and long rise favors higher mass. If we assume a smooth rise, the best-fit ejecta mass is 15 M${}_{\odot }$. The slow tail phase is well-matched due to a long spin-down time, given by the relatively weak B-field (Kasen & Bildsten 2010; Nicholl et al. 2013). This model gives a fairly good match to the velocity evolution, especially the slow velocity decline at late times. On the other hand, the broad luminosity peak in this case is a rather poor fit to the data.

Zoom In Zoom Out Reset image size

Table 1.  Magnetar Fit Parameters for Figure 19

Rise time/day	${M}_{{\rm{ej}}}$/M${}_{\odot }$ a	P/ms	B/1014 G	${\chi }_{{\rm{red}}}^{2}$
69b	8.4	2.1	0.9	0.31
85	15.1	1.7	1.0	0.48

Notes.

^aAssuming $\kappa =0.2$ cm² g⁻¹. ^bEstimated from scaled light curve of LSQ14bdq.

Download table as:  ASCIITypeset image

If we instead assume that the first point is part of the bump phase, we get a more convincing fit between −50 days and maximum, with a lower ejecta mass of ${M}_{{\rm{ej}}}$ ≈ 8.4 M${}_{\odot }$. The model still peaks a little earlier than the data, which is interesting in light of the comparison to PTF12dam in the lower panel of Figure 18. This model also gives an excellent match to the photospheric velocity. While the temperature in this model fits the data well from −30 to +50 days, it is too low by up to a factor of ∼2 by 250 days. However, the fact that the better overall fit to the light curve, temperature and velocity is the one for which we assume an initial bump phase is consistent with the analysis in Section 3.3.

However, inspection of Table 2 shows that this model has a total rise time of 69 days, seemingly in contradiction to the presumed bump detection at −79 days. The widths of the bumps may be 10–15 days (Nicholl et al. 2015a; Nicholl & Smartt 2016), which could then imply a discrepancy between the explosion time and the magnetar fit of up to 20–25 days. There are two possibilities. The first is that this is a real inconsistency, and the magnetar fit must be fixed with an explosion time at the first detection at the latest (i.e., our "slow rising" model). The second is that thermalization of the magnetar wind energy in the ejecta is not 100% efficient from the time of explosion. Our model does naively assume that the magnetar energy input is 100% efficient but this is by no means certain. The detailed physical processes by which the energy is emitted and thermalised are not yet well understood (Kasen et al. 2015; Metzger et al. 2015). Therefore there could be a delay between core-collapse and the onset of efficient energy transfer between the magnetar and the expanding ejecta. Given these uncertainties, we conclude that a magnetar spin-down model with relatively high ejecta mass and low B-field is a good candidate to explain the properties of SN 2015bn.

Table 2.  Interaction Fit Parameters for Figure 20

Rise time/day	Profile	${M}_{{\rm{ej}}}$/M${}_{\odot }$a	${M}_{{\rm{CSM}}}$/M${}_{\odot }$ a	${M}_{{\rm{Ni}}}$/M${}_{\odot }$	${E}_{{\rm{k}}}/{10}^{51}$ erg	${R}_{{\rm{int}}}$/1015cmb	$\mathrm{log}({\rho }_{0}$/g cm−3)c	${\chi }_{{\rm{red}}}^{2}$
67	shelld	49.0	18.5	⋯	1.9	3.0	−10.94	0.51
87	shell	21.9	19.2	⋯	1.8	3.0	−12.00	0.35
82	shell	13.8	13.0	3.7	1.3	3.0	−12.34	0.15
61	winde	30.1	15.1	⋯	4.3	0.3	−10.49	0.37
85	wind	33.3	18.0	⋯	4.5	0.3	−10.61	0.40
59	wind	29.4	16.1	2.3	4.4	3.1	−12.19	0.17

Notes.

^aAssuming $\kappa =0.2$ cm² g⁻¹. ^bInteraction radius = inner radius of CSM density profile. ^c ${\rho }_{0}={\rho }_{{\rm{CSM}}}(r={R}_{{\rm{int}}})$. ^d ${\rho }_{{\rm{CSM}}}(r)={\rho }_{0}$. ^e ${\rho }_{{\rm{CSM}}}(r)={\rho }_{0}{r}^{-2}$.

Download table as:  ASCIITypeset image

Next we investigate models powered by a collision between the fast SN ejecta and dense CSM. In this context, the observed luminosity comes from reprocessing and thermalization of the kinetic energy by a forward shock wave propagating in the CSM and a reverse shock in the ejecta. For the energy conversion to be efficient, in general the mass contained in the CSM must be an appreciable fraction of that ejected by the explosion—i.e., at least a few solar masses (Chatzopoulos et al. 2013; Nicholl et al. 2014; Inserra et al. 2015). The progenitors of interacting SLSNe would thus require exceptionally high mass-loss rates, perhaps as much as 0.1–1 M${}_{\odot }$ yr⁻¹ (e.g., Benetti et al. 2014). This is difficult to achieve via stellar winds, and could instead point toward large discrete mass ejections, such as the eruptions of luminous blue variables (Justham et al. 2014) or hypothetical pulsational-PISNe (Woosley et al. 2007). Thus the expected CSM density profile around the star is quite uncertain. Following Chatzopoulos et al. (2012), we calculate models with two representative density profiles: ${\rho }_{{\rm{CSM}}}\propto {r}^{-2}$, as expected for a wind; and ${\rho }_{{\rm{CSM}}}=$ constant, which may be more appropriate for a CSM "shell" from a discrete mass ejection. Because the interaction model has significantly more free parameters than the magnetar models, there is no tight relationship between the rise and decline rates (Nicholl et al. 2015b). Therefore it is not so important in this case how we treat the first detection of SN 2015bn at −79 days; typically, increasing the ejecta or CSM mass can extend the rise time sufficiently to include this point as desired, without significantly altering the shape during the decline phase. For each density profile, we calculate a slow- and fast-rising model to demonstrate this (see fit parameters in Table 2).

The fits are shown in Figure 20. In the top panel we plot models with uniform shells of CSM. We fit for ${M}_{{\rm{ej}}}$, ${M}_{{\rm{CSM}}}$, ${\rho }_{{\rm{CSM}}}$, the radius of the ejecta-CSM interface and kinetic energy (${E}_{{\rm{k}}}$). We note that the radius, ${R}_{{\rm{int}}}$, is only weakly constrained by the light curve fitting routine, but can be checked against the observed temperature evolution³⁶ (bottom panel). The best-fit slow-and fast rising shell models both have ${M}_{{\rm{CSM}}}$ ≈ 19 M${}_{\odot }$ and ${E}_{{\rm{k}}}$ ≈ 2 × 10⁵¹ erg. The ratio ${M}_{{\rm{ej}}}$/${M}_{{\rm{CSM}}}$ is of order unity for the model with the slow rise, while the faster rising model has ${M}_{{\rm{ej}}}$/${M}_{{\rm{CSM}}}$ = 2.6. These are fairly similar to the ratios found by Nicholl et al. (2014) for a number of other SLSNe I. However, these curves are unable to reproduce the transition to the shallower decline rate after 60 days. By maximum light, the forward and reverse shocks have finished traversing the CSM and ejecta, respectively, and therefore no further energy can be deposited. The light curves therefore decline exponentially from peak as the stored energy diffuses out.

Zoom In Zoom Out Reset image size

The middle panel shows models using the wind density profile. In general, the wind scenario favors similar or slightly lower ${M}_{{\rm{ej}}}$ and ${M}_{{\rm{CSM}}}$ compared to the shell models, with ${M}_{{\rm{ej}}}$ ≈30 M${}_{\odot }$ giving a good fit along with ${M}_{{\rm{CSM}}}$ = 15–18 M${}_{\odot }$ to cover the range of possible rise times. The kinetic energy in the wind model is 4 × 10⁵¹ erg. The wind models are more sensitive to the interaction radius, ${R}_{{\rm{int}}}$, compared to the shell models. For the inferred radii of $\approx 3\times {10}^{14}$ cm, a wind speed of 10 km s⁻¹ would imply that the mass loss occurred over a period of 10–200 days before explosion. However, these models suffer from the same problem at late times as do the shell models, in that they decline too quickly. Comparing the temperature evolution of each model (bottom panel), we find that both the wind and shell models give a reasonable match to observations, with the slow rising models predicting cooler temperatures that are closer to the data. While the quality of the light curve fits do not allow us to distinguish between wind and shell models, there exists significantly more literature describing the interaction of SN ejecta with a wind density profile. Therefore the remaining analysis will focus more on the wind model. This will also be useful when we come to model the radio light curve in Section 7.

No explicit velocity comparison is possible, since the simple interaction model assumes a stationary CSM. However, as is often the case when modeling SLSNe, the exclusively high-velocity lines seen in the spectrum may be difficult to reconcile with having $\gtrsim 10$ M${}_{\odot }$ of slower-moving material around the star. In an attempt to address this issue more quantitatively, we can use some of the analytic equations given by Chevalier & Irwin (2011) in conjunction with our fits. Those authors parameterised the wind as a function of radius as ${\rho }_{w}\;=\dot{M}/(4\pi {v}_{w}){r}^{-2}\equiv 5.0\quad \times \quad {10}^{16}{D}_{* }{r}^{-2}$, where ${\rho }_{w}$ is the density of the wind, $\dot{M}$ is the mass-loss rate, v_w is the wind velocity, and

$$\begin{eqnarray}&&{D}_{* }\equiv \displaystyle \frac{\dot{M}}{{10}^{-2}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}{\left(\displaystyle \frac{{v}_{{\rm{wind}}}}{10\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-1}.\end{eqnarray} \tag{ 1 }$$

We found that the peak luminosity is quite sensitive to ${D}_{* };$ the value for the wind models in Figure 20 is ${D}_{* }=19.4$ (or equivalently, $\dot{M}/{v}_{10}\approx 0.2$ M${}_{\odot }$ yr⁻¹, where v₁₀ is the wind velocity in units of 10 km s⁻¹). In the Chevalier & Irwin (2011) framework, this can be used to calculate an effective diffusion radius:

$$\begin{eqnarray}&&{R}_{d}=4.0\times {10}^{14}\;{\kappa }_{0.34}^{0.8}\;{E}_{{\rm{k,51}}}^{0.4}\;{M}_{{\rm{ej,10}}}^{-0.2}\;{D}_{* }^{0.6},\end{eqnarray} \tag{ 2 }$$

where the subscripts refer to normalization of the parameters to 0.34 cm² g${}^{-1},{10}^{51}$ erg and 10 M${}_{\odot }$ as used by those authors. Our light curve models give estimates for all of these parameters (Table 2), leading to ${R}_{d}=7.3\times {10}^{15}$ cm, which is reassuringly close to the radii determined from blackbody fits (Figure 17).

Chevalier & Irwin (2011) found that the SN behavior depended on whether the terminal radius of the wind, R_w, was larger or smaller than R_d. They predicted that for ${R}_{w}\gt {R}_{d}$, a dense shell forms deep in the wind, leading to low velocities, and continued late-time interaction that gives a flattening in the light curve. On the other hand, if ${R}_{w}\lt {R}_{d}$, the outermost layers are accelerated by radiation pressure, and there is little power input after maximum light. The authors suggest that SLSNe I are examples of the latter situation, whereas SLSNe IIn correspond to the former. For SN 2015bn, our models give R_w = 5.6–6.4 × 10¹⁵ cm, depending on the CSM mass. Thus we are in the regime ${R}_{w}\lesssim {R}_{d}$, possibly resolving the lack of narrow lines in the spectrum. In this case we do not expect to see a flattening of the light curve, which is indeed borne out by our model fits (though not by the data). Therefore another energy source must power the tail. ${R}_{w}\lesssim {R}_{d}$ was also found to be the case for SN 2010gx (Chevalier & Irwin 2011); however Chomiuk et al. (2011) applied the same analysis to PS1-10ky and PS1-10awh and found ${R}_{w}\gtrsim {R}_{d}$. Within this approximate framework, the more robust constraint seems to be that ${R}_{w}\approx {R}_{d}$, which is required for efficient production of radiation (Chomiuk et al. 2011). Finally, the Chevalier & Irwin (2011) model allows us to calculate a shock velocity of $\approx 7000$ km s⁻¹ assuming an 80 day rise time based on the observed data.

In order to recover the slow decline observed after +50 days, we re-ran the fits with ${M}_{{\rm{Ni}}}$ as an additional free parameter, and found that adding 2–4 M${}_{\odot }$ of ⁵⁶Ni to each model could reproduce the late-time behavior. This is shown in Figure 20. For the wind model, a lower ${\rho }_{{\rm{CSM}}}$ then gave a match the light curve shape around peak, with a corresponding increase in the interaction radius to keep the peak brightness constant (this is equivalent to keeping the pre-explosion mass-loss rate and velocity constant). The interaction models with ⁵⁶Ni-decay give good fits to the light curve, although the values of ${M}_{{\rm{Ni}}}$ are extremely large—normal core-collapse SNe have ${M}_{{\rm{Ni}}}$ ≈ 0.1 M${}_{\odot }$, and even LGRB-SNe synthesize ${M}_{{\rm{Ni}}}$ ≲ 0.5 M${}_{\odot }$ (e.g., see the recent review by Cano et al. 2016).

The alternative scenario sees the shallower tail phase of the light curve powered by further interaction with mass-loss at a larger radius. This is appealing, since we do not need to invoke both a large CSM mass and a large ⁵⁶Ni mass. We construct such a model by first fitting an interaction model to the data at $t\lt 60\;{\rm{days}}$, and then subtracting this from the whole light curve. We then fit the residuals at $t\gt 60\;{\rm{days}}$ as an interaction with additional CSM. Most of the model parameters for this second interaction are fixed by the fit parameters of the first interaction. In the following discussion, we use the subscripts 1 and 2 to refer to the fits around peak and to the late-time residuals, respectively. The fixed parameters are ${M}_{{\rm{ej,2}}}={M}_{{\rm{ej,1}}}+{M}_{{\rm{CSM,1}}};$ ${E}_{{\rm{k,2}}}={E}_{{\rm{k,1}}}\;-{E}_{{\rm{rad,1}}}$ (where ${E}_{{\rm{rad}}}$ is the total energy lost through radiation, $\approx 2\times {10}^{51}$ erg); ${R}_{{\rm{int,2}}}={R}_{{\rm{CSM,1}}}+\sqrt{(10{E}_{{\rm{k,2}}})/(3{M}_{{\rm{ej,2}}})}{\rm{\Delta }}t$ (where ${R}_{{\rm{CSM}}}$ is the outer radius of the CSM, and ${\rm{\Delta }}t$ is the time between shock breakout from the inner CSM and the second interaction—this is approximately equal to the time from maximum light until the beginning of the tail phase, i.e., $\approx 60$ days). Therefore the only free parameters are ${M}_{{\rm{CSM,2}}}$ and ${\rho }_{\mathrm{0,2}}$. For a given density, the CSM mass only affects the duration of the second interaction but not the luminosity, thus this parameter is only weakly constrained by the data. However, we find that we require ${M}_{{\rm{CSM,2}}}\gtrsim 10$ M${}_{\odot }$.

The best-fitting model is shown in Figure 21. While the late-time fit is not as good as the ⁵⁶Ni-powered model in Figure 20, we caution that we are now pushing the simplified analytic CSM model to its limits, and this fit is intended only to give an order-of-magnitude estimate for the mass-loss rate required to give the shallower tail phase. In this case we find $\dot{M}/{v}_{10}\approx 0.05$ M${}_{\odot }$ yr⁻¹. (See Table 3 for a full description of the model parameters). As this model is only approximate, it is possible that a continuous CSM could also fit the data, for example a dense inner shell, attached to an outer wind component. Regardless, the luminosity at late times seems to require a mass-loss rate $\dot{M}/{v}_{10}\gtrsim 0.01$ M${}_{\odot }$ yr⁻¹. However, this model raises questions about the spectroscopic evolution. While it is conceivable that the single-shell model could avoid showing narrow lines in certain circumstances, it seems likely that the outer CSM, which remains unshocked around light curve maximum, would imprint low-velocity features on the observed spectrum. This is essentially the same conclusion reached by Chevalier & Irwin (2011), as discussed above. One possible way to avoid seeing narrow lines at early times could be if this outer CSM were a face-on disk or torus.

Zoom In Zoom Out Reset image size

Table 3.  Two-component Interaction Fit Parameters for Figure 21

	Profile	${M}_{{\rm{ej}}}$/M${}_{\odot }$ a	${M}_{{\rm{CSM}}}$/M${}_{\odot }$ a	${M}_{{\rm{Ni}}}$/M${}_{\odot }$	${E}_{{\rm{k}}}/{10}^{51}$ erg	${R}_{{\rm{int}}}$/1013 cmb	$\mathrm{log}({\rho }_{0}$/g cm−3)c	${\chi }_{{\rm{red}}}^{2}$
1	windd	25.0	15.0	⋯	4.0	140	−12.3	1.13
2	wind	40.0e	15.0	⋯	2.0f	350	−14.1	0.26

Notes.

^aAssuming $\kappa =0.2$ cm² g⁻¹. ^bInteraction radius = inner radius of CSM density profile. ^c ${\rho }_{0}={\rho }_{{\rm{CSM}}}(r={R}_{{\rm{int}}})$. ^d ${\rho }_{{\rm{CSM}}}(r)={\rho }_{0}{r}^{-2}$. ^e ${M}_{{\rm{ej,2}}}={M}_{{\rm{ej,1}}}+{M}_{{\rm{CSM,1}}}$. ^f ${E}_{{\rm{k,2}}}={E}_{{\rm{k,1}}}-{E}_{{\rm{rad,1}}}$.

Download table as:  ASCIITypeset image

Table 8.  Publicly Available CRTS Photometry

MJD	Phasea	Rb
57014.00	−79.0	18.87 (0.28)
57064.00	−34.1	16.95 (0.18)
57091.44	−9.5	16.66 (0.04)
57097.44	−4.1	16.63 (0.02)
57102.44	0.4	16.65 (0.02)
57109.44	6.7	16.65 (0.04)
57123.44	19.3	16.88 (0.02)
57129.44	24.6	17.04 (0.02)
57135.44	30.0	17.15 (0.05)
57155.44	48.0	17.33 (0.12)
57162.44	54.3	17.44 (0.11)

Notes.

^aIn rest-frame days from MJD = 56102. ^bAverage magnitude of $\approx 4$ detections on each night.

Download table as:  ASCIITypeset image

In Section 5, it was observed that the late-time decline rates of SN 2015bn and a number of other SLSNe I look very similar to the radioactive decay of ⁵⁶Co, the daughter nucleus of ⁵⁶Ni. We also saw above that introducing a few solar masses of ⁵⁶Ni to the interaction models gave a good fit to the tail luminosity. However, powering this phase by ⁵⁶Co alone would imply that $\approx 15$ M${}_{\odot }$ of ⁵⁶Ni were synthesized in the explosion, using the scaling relations of Hamuy (2003). This would require the pair-instability explosion of a stellar core of over 110 M${}_{\odot }$, according to the models of Kasen et al. (2011). Such a scenario was ruled out for PTF12dam by Nicholl et al. (2013), based on the rise time and blue spectra. In addition to the spectroscopic and photometric similarity between these two events, the undulations in the light curve of SN 2015bn seem to be incompatible with a PISN, as such a huge ejected mass would be expected to wash out any undulating structure in the light curve due to the very long diffusion time (e.g., Fraser et al. 2013). Thus it is difficult to conceive of a plausible model in which the tail phase of SN 2015bn is predominantly powered by ⁵⁶Co-decay. However, the conspiracy among an increasing number of SLSNe I to match this decline rate is certainly intriguing, and we therefore examine radioactively powered models in some detail here.

We fitted the light curve of SN 2015bn with a ⁵⁶Ni-powered model, finding that reasonable fits were only possible with ${M}_{{\rm{Ni}}}$ ≈ ${M}_{{\rm{ej}}}$. The best-fit model shown in Figure 22 has ${M}_{{\rm{ej}}}$ ≈ 27.0 M${}_{\odot }$, ${M}_{{\rm{Ni}}}$ ≈ 26 M${}_{\odot }$. Note that only one fit is shown, as the best-fitting models including/excluding the first data point are almost identical. In order to decrease the diffusion time to fit the initial decline from peak, a large explosion energy of ${E}_{{\rm{k}}}$ = 3 × 10⁵² erg was required. This in turn reduced the γ-ray trapping efficiency at late times, leading to a final decline rate that is faster than ⁵⁶Co-decay. The model is unconvincing both in terms of the unrealistic ⁵⁶Ni-fraction in the ejecta and a poor overall fit to the light curve. Moreover, this model predicts a low photospheric temperature, and a very high velocity. The high velocity is due to the large value of ${E}_{{\rm{k}}}$. In general, this fit is inferior to the magnetar and interaction models.

Zoom In Zoom Out Reset image size

We also compare to hydrodynamic simulations of PISNe calculated by Kasen et al. (2011). The peak luminosity is well-matched by the explosion of a 130 M${}_{\odot }$ stripped-envelope PISN, synthesising 40 M${}_{\odot }$ of ⁵⁶Ni. It also matches the decline rate at times $\gtrsim 50$–100 days, but the faster decline observed in the data from 0 to 50 days is difficult to reconcile with this model. The rise time of this model is also too long relative to the data, even including the earliest point. The absolute luminosity of SN 2015bn at this phase is in better agreement with the 120 M${}_{\odot }$ model, which has 24 M${}_{\odot }$ of ⁵⁶Ni, similar to our own diffusion fit (the decay energy in PISN models remains fully trapped at late times due to the large ejected mass). However, this lower-mass model severely underestimates the peak luminosity and does not reproduce the behavior for times earlier than +50 days. One possibility that could perhaps match the light curve shape would be to add strong CSM interaction around peak to the 120 M${}_{\odot }$ PISN model. However, this would require invoking an extreme CSM mass as well as the extreme ejecta and ⁵⁶Ni masses, and a detailed investigation of such a model is beyond the scope of this paper. The temperature in the more massive (hotter) 130 M${}_{\odot }$ model is cooler than SN 2015bn by about 3000 K at all times. The velocity in this model gives a fairly good match until about +50 days, but the rapid decline as the photosphere recedes looks quite different from the slow decline exhibited by SN 2015bn.

In Figure 23, we compare the +243 day spectrum of SN 2015bn to the model spectrum of a 130 M${}_{\odot }$ hydrogen-poor PISN, at a similar phase, from Jerkstrand et al. (2016). This model differs from the nebular PISN model shown by Gal-Yam et al. (2009), because the Jerkstrand et al. model is computed using realistic ejecta models from Heger & Woosley (2002), whereas the Gal-Yam et al. models used tuned parameters to fit the spectrum of SN 2007bi, resulting in relative abundances that did not match PISN explosion models (Nicholl et al. 2013). The new models also include line blocking, which is extensive in these high-density ejecta and damps UV/blue emission. The model shown gives a good match to SN 2015bn between ∼6000–8000 Å (oxygen and calcium). We note that strong oxygen and calcium lines in the nebular phase are not unique to PISNe, but are generically seen in core-collapse SNe (e.g., Filippenko 1997). The match to the line widths and overall luminosity scale are good. This match at red wavelengths is much better than was seen in similar comparisons between the PISN models and SN 2007bi and PTF12dam (Jerkstrand et al. 2016). However, there is a huge discrepancy between the model and data below $\approx 5500\;\mathring{\rm A}$. In the PISN model, line-blanketing by the large iron mass suppresses the flux in the blue, whereas the late spectrum of SN 2015bn retains a pronounced blue continuum that is brighter than the model by an order of magnitude. This discrepancy has previously been pointed out for other slow-declining SLSNe by Dessart et al. (2012), Nicholl et al. (2013), Jerkstrand et al. (2016). We note that line formation at bluer wavelengths may be more complex compared to in the red, but it is unclear what model deficiencies could account for the scale of this difference between model and data. Despite the similarity in the light curve slope and spectrum above 6000 Å, there remain many important differences between the observational data and models powered exclusively by ⁵⁶Ni.

Zoom In Zoom Out Reset image size

Having discussed the possible power sources for the bulk of the light curve, we now turn to interpretations of the shoulder and knee. In both cases, the temperature evolution suggests that an additional heat source is active over these times, as all of the power sources we have considered (⁵⁶Ni-decay, magnetar, or shock heating from interaction) predict a monotonic decline from maximum. The discussion in this section will focus on the two models that gave the more convincing fits in the previous section (magnetar spin-down and CSM interaction), beginning with interaction.

Similar light curve undulations to those observed in SN 2015bn were also seen in the light curve of SN 2009ip, which is thought to be either a SN IIn or a SN imposter (Pastorello et al. 2013). Either way, the luminosity at maximum light was powered by interaction with pre-expelled material (e.g., Fraser et al. 2013; Mauerhan et al. 2013; Smith et al. 2013; Graham et al. 2014; Margutti et al. 2014a). However, Fraser et al. (2013) note that an interaction between more than a few solar masses of ejecta and CSM (which would certainly be required to match the high luminosity in SN 2015bn) would struggle to match the undulations, as these would tend to be washed out by the long diffusion time. Therefore the most probable way to explain these using CSM interaction is to introduce multiple collisions with additional shells or clumps of material.

All of the aforementioned studies of SN 2009ip found that at least some of the features in the light curve could be explained by successive collisions with mass expelled in previous observed outbursts of the progenitor. If we believe that the historical CSS detection of SN 2015bn at $\approx 185\;{\rm{days}}$ before explosion was an outburst of the progenitor, and that interaction with this material powers the knee or shoulder, the blackbody radii at these epochs (Figure 17) would imply that the pre-expelled material had a velocity of ∼3000 km s⁻¹. This would be compatible with pulsational pair-instability mass ejections (Woosley et al. 2007).

However, the blackbody radius in SN 2015bn immediately begins to decrease at the onset of the knee phase. This may be difficult to explain in an interaction model, as the collision with a massive shell should initially be optically thick (e.g., Moriya & Maeda 2012). In SNe powered by interaction, the radius of the emitting surface generally shows an increase as the shock propagates through the CSM (e.g., Fassia et al. 2000; Margutti et al. 2014a). This may not be too problematic, as only a relatively small amount of CSM is needed to power the undulations. For example taking the knee feature, we find that SN 2015bn is over-luminous by 2 × 10⁴³ erg s⁻¹ for 15 days, compared to a smooth decline. Assuming a velocity v = 7000 km s⁻¹ from the spectrum at this phase, we can use the common scaling relation (Quimby et al. 2007; Smith & McCray 2007) $L\approx \tfrac{1}{2}{M}_{{\rm{CSM}}}{v}^{2}/{t}_{{\rm{rise}}}$, giving ${M}_{{\rm{CSM,knee}}}\sim 0.05$ M${}_{\odot }$. Such a low CSM mass is probably not compatible with pulsational-PISN shells; the final pulses before explosion in all of the Woosley et al. (2007) models eject at least $\gtrsim 1$ M${}_{\odot }$. The low mass, relative to the ejecta/CSM masses in any of our model fits, might help to explain why we see little change in the continuum (with the temperature approximately constant at this phase), but this modest CSM mass may be expected to result in narrow emission lines, as the optical depth will be low compared to the case of massive CSM. However, as for other SLSNe I, no narrow lines are observed in the spectrum of SN 2015bn.

Next we consider the magnetar model. Before investigating any magnetar-specific means of producing the undulations in the light curve, it should be pointed out that the low-mass CSM interaction described above could equally apply to a light curve that was primarily powered by magnetar spin-down. Nevertheless, the magnetar scenario does have other means of producing fluctuations in the light curve. One such mechanism is the magnetar-powered shock breakout described by Kasen et al. (2015). In this model, the central overpressure from the magnetar wind drives a second shock through pre-exploded ejecta, which breaks out at large radius and hence can give a bright optical display. However, their Equation (26) shows that this breakout should happen within 20 days of explosion for any sensible combination of ejecta mass, explosion energy and magnetar parameters. Moreover, the breakout would not be noticeable in the light curve if it occurred near maximum brightness. Therefore this model is only applicable to early-time bumps, and not to the shoulder or knee in SN 2015bn.

Metzger et al. (2014) found that for a magnetar-powered SN, the hard radiation field should drive ionization fronts outwards from the base of the ejecta. Being the most abundant element in SN Ic ejecta, oxygen was considered to be the dominant source of electrons at early times. Their model predicted that when the O ii layer reached the ejecta surface, the opacity to UV photons would decrease, leading to a UV-breakout. This occurred on a timescale of tens of days after explosion, and this timescale increases with ejecta mass. Comparing with the properties of SN 2015bn during the shoulder phase, this model seems to consistently explain its timing, the fact that it is only apparent in the blue and UV bands, and the presence of O ii lines in the spectrum at this time. The fast decline in the color temperature during this phase also indicates that it could be a form of breakout event.

A key prediction of this model is that the UV breakout should be followed by X-ray breakout tens to hundreds of days later. This occurs when the O iv ionization front reaches the surface. For the model computed by Metzger et al. with a 2 ms spin period (needed to power the peak optical luminosity of SN 2015bn), X-ray breakout occurred 240 days after explosion. However, that model assumed an ejecta mass of only 3 M${}_{\odot }$, a factor of 3–5 lower than the masses inferred from our fits to SN 2015bn (Table 1). The breakout timescale is a strong function of ejecta mass, and thus X-ray breakout is not expected to occur over the timescale of our observations. However, we are carrying out an X-ray monitoring campaign for SN 2015bn to look for signatures of ionization breakout; the results will be published in a future study by Margutti et al.

During the knee plateau from +30 to 50 days, the photospheric radius starts to decrease, and equivalent widths of spectral lines from heavier species such as iron and silicon increase. This seems to indicate that the photosphere is beginning to undergo significant recession into the ejecta. The velocity evolution slows at around this time, suggesting that the velocity profile is relatively flat within a certain radius. Kasen & Bildsten (2010) predicted that the magnetar wind would sweep up most of the ejecta into a dense shell with uniform velocity. Moreover, they show that the temperature jumps sharply at the edge of this shell. The knee, and the associated increase in temperature and flattening in velocity, could therefore be explained by the photosphere crossing into this hotter region as it recedes. If SLSNe I are generally powered by magnetars, one may wonder why this distinctive behavior is not seen more often. The answer could be that, in lower mass events, the magnetar wind sweeps up essentially all of the ejecta, such that there is no discontinuity. Alternatively, it could simply be that most SLSN light curves do not have sufficient temporal sampling for us to catch variations on this timescale.

Another possibility to be considered is that one or both of the undulations are powered by recombination, similar to the plateaus in SNe IIP. However, the high temperature and the spectra clearly rule out hydrogen recombination. One plausible candidate that could be consistent with the shoulder is oxygen recombination. Oxygen in SNe is expected to be mostly singly ionized at temperatures between 12,000 and 15,000 K, and neutral below this (Hatano et al. 1999; Quimby et al. 2013). SN 2015bn is seen to drop below ∼12,000 K at precisely this phase. The early spectroscopic evolution may also support this: the O ii lines in SN 2015bn are weak at ∼20 days pre-peak compared with those in PTF12dam and other SLSNe at similar phases (see Figures 13 and 14). Instead, the early spectrum closely matches PTF12dam at around maximum light, when these lines are close to disappearing. This could suggest that O ii is being depleted by recombination. The photospheric recession during the knee phase may also point toward recombination. However, it is unclear in this scenario why the features would be most pronounced in the UV bands, or which ion could recombine to power the knee (${T}_{{\rm{BB}}}\approx 8500$ K).

Undulations or plateaus in the light curve evolution could also in principle result from an abrupt change in the continuum optical/UV opacity, as might arise due to the sudden emergence of an ionization front through the ejecta (Metzger et al. 2014). If the ejecta remain largely neutral at early times, the optical opacity will be due primarily to bound-bound transitions and could be relatively modest ($\kappa \lesssim 0.05$ cm² g⁻¹, depending on the composition). However, once the ejecta become ionized and the number of free electrons increases, the resulting electron scattering opacity would come to dominate, increasing κ by a factor of several. This sudden rise in opacity could slightly delay the escape of radiation from the magnetar nebula or CSM interaction shocks, producing a plateau or lull in the light curve decline. Exploring this scenario further would require a detailed radiative transfer calculation, which is beyond the scope of this paper.

To see the undulations more clearly, and to look for any hints of periodicity, we removed the gross structure of the light curve to examine the residuals. We did this by fitting and then subtracting low-order polynomials (we tried first-, second- and third-order). We also experimented with subtracting one of the magnetar-powered fits (solid line in Figure 19). Figure 24 shows the bolometric light curve after subtraction of the various fits. The detrended data are seen to exhibit variation on a timescale of $\approx 30$–50 days. A dominant timescale of $\lesssim 50\;{\rm{days}}$ would be inconsistent with a very massive model, such as a PISN. However, it is fully consistent with the diffusion timescale in an ejecta of $\approx 7$–15 M${}_{\odot }$ (Arnett 1982; Nicholl et al. 2015b), which is similar to what we infer from our magnetar fits. Thus the undulations seem to be compatible with variations in the velocity/density structure in ejecta similar to that in our magnetar models.

Zoom In Zoom Out Reset image size

If the undulations are really periodic, this could indicate that they are caused by interaction with a spiral density perturbation, perhaps caused by binary interaction (Fraser et al. 2013). Evidence of close binary interaction would certainly be an intriguing part of the picture for understanding SLSN progenitors. Unfortunately, the data do not span a sufficiently long time interval (Δt ≈ 150 days) to reliably pick out periodicity on a timescale on this order (t₀ ≈ 50 days ≈ Δt/3). For example, Martin et al. (2014) applied a periodogram analysis to SN 2009ip, but ignored peaks in the power spectrum with periods $\gtrsim {\rm{\Delta }}t/3$. Although the undulations in the light curve of SN 2015bn are an important clue to the nature of this explosion, sampling over much longer timescales will be needed to robustly test for periodicity or the presence of a dominant timescale.

We first tested whether these limits can constrain the presence of an off-axis γ-ray burst associated with SN 2015bn. Lunnan et al. (2014) found that SLSNe I and LGRBs tend to occur in similar environments, though recently Angus et al. (2016) have proposed that the only similarity is low metallicity. In any case, both types of explosions display similar spectra lacking in hydrogen and helium (e.g., Pastorello et al. 2010), and both may be powered by central engines. Recently, Greiner et al. (2015) presented a luminous ($M\approx -20$) SN associated with an ultra-long GRB, but it was not clear from the relatively featureless spectrum whether this was a "classical" SLSN I. The presence of a radio afterglow from a confirmed SLSN would firmly establish the SLSN-GRB connection; alternatively, a lack of radio emission could instead indicate that either no stable jet forms, or it does not break out of the massive progenitor, and thus does not accelerate the outer ejecta to relativistic velocity.

We compare models for LGRB-SNe at a variety of viewing angles in Figure 25. The simulations are described by Van Eerten et al. (2010), and assume typical LGRB parameters: total kinetic energy in the two jets of ${E}_{{\rm{k,jets}}}=2\times {10}^{51}$ erg; a homogeneous particle density of n = 1 cm⁻³ in the surrounding medium; equal fractions of the total energy density in magnetic fields and relativistic electrons (${\epsilon }_{B}={\epsilon }_{e}=0.1$); jet half-opening angle ${\theta }_{{\rm{jet}}}=0.2\;{\rm{rad}};$ isotropic-equivalent energy ${E}_{{\rm{iso}}}={10}^{53}$ erg; and a power-law slope of p = 2.5 for the accelerated electron distribution. For explosions viewed further from the jet axis, the radio emission is weaker, and peaks later. The models are calculated at a frequency of 8.5 GHz, which gives a good match to the observer-frame 7.4 GHz limit for SN 2015bn.

Zoom In Zoom Out Reset image size

Our deep radio limit for SN 2015bn rules out even the most off-axis model, making SN 2015bn the first SLSN for which we can robustly exclude a luminous LGRB. However, we also compare to a sample of observed "low-luminosity" (ll) GRBs, for which the emission does not seem to be strongly collimated as in LGRBs. ll-GRBs may form a separate population from high-luminosity LGRBs, and dominate the volumetric GRB rate locally (Guetta & Della Valle 2007; Liang et al. 2007). As demonstrated by the figure, we cannot rule out a ll-GRB associated with SN 2015bn, and can only exclude the luminous and highly collimated jets of a "healthy" LGRB. Similarly deep radio observations of the next nearby SLSN I before or around optical maximum should be able to constrain the presence of a ll-GRB.

Radio observations can also be useful in breaking the degeneracies between the competing models for the optical light curve. While the magnetar and interaction models give similarly good fits to the optical data, their predicted radio signatures are very different. In general, radio synchrotron emission is generated by electrons accelerated to relativistic velocities at the forward shock (Chevalier 1982; Weiler et al. 2002; Chevalier & Fransson 2006).

If SLSNe I are otherwise fairly normal SNe Ic that are reheated by a magnetar engine, we expect the radio properties to resemble those of typical SNe Ic. For those objects, the CSM density is relatively low, with radio light curves implying mass-loss rates of around 10⁻⁷–10⁻⁵ M${}_{\odot }$ yr⁻¹ (if the wind velocity is 10³ km s⁻¹; e.g., Berger et al. 2002; Soderberg 2007; Drout et al. 2015). The interaction is dominated by the outermost, fastest ejecta ($v\gtrsim 0.15c$), which carries only ∼10⁻⁵ of the total kinetic energy (Matzner & McKee 1999), while the dominant source of absorption is synchrotron self-absorption (Chevalier 1998). Most normal SNe Ic (i.e., no LGRB) that have been detected in the radio peaked at ${L}_{\text{8 GHz}}\lesssim {10}^{28}$ erg s⁻¹ Hz⁻¹, on a timescale of tens of days (Berger et al. 2003; Soderberg et al. 2010; Margutti et al. 2014b), though there are some exceptions (e.g., Corsi et al. 2014; Kamble et al. 2016). Generally, the predicted luminosity at ∼300 days is well below our limit for SN 2015bn. Radio emission from the magnetar itself is not expected, as Galactic magnetars are not detected as radio sources (e.g., Gaensler et al. 2001), and any signal that was emitted would be obscured by free–free absorption in the ionised inner regions of the ejecta (Chomiuk et al. 2011). In conclusion, our radio observation is not constraining for the magnetar model.

For the alternative model in which the optical luminosity is powered by interaction with CSM, the masses and densities involved are much higher; our best-fit wind model has a mass-loss rate higher than that inferred for normal SNe Ic by several orders of magnitude (if normalized to the same wind velocity). In this regime, all of the ejecta interacts with the massive CSM, and therefore the kinetic energy involved is the total energy of the ejecta, $\approx {\rm{few}}\times {10}^{51}$ erg (Section 6). For the very dense mass-loss in our CSM fits, we expect that the radio emission will be obscured by free–free absorption until well after maximum light (for example radio emission in SNe II—which have mass-loss rates intermediate between SNe Ic and our models here—peaks at $\gtrsim 100\;{\rm{days}}$ after explosion; Weiler et al. 2002).

For each of the models shown in Figure 20, the forward shock reaches the outer edge of the CSM at the time of optical maximum; the declining light curve is then powered by diffusion of the stored shock energy. This is similar to the statement that ${R}_{w}\approx {R}_{d}$ (Section 6.3). Nicholl et al. (2014) fitted a number of SLSN light curves with this interaction model, and also found that good fits required shock breakout at around maximum light. In this scenario, we might expect to see radio emission around optical maximum, when the shock is near the outer surface of the CSM. Once the shock expands past the CSM, there is no further particle acceleration. Adiabatic expansion rapidly reduces the energy density in relativistic particles and magnetic fields, and the radio luminosity falls off rapidly (${L}_{\nu }\propto {t}^{-6};$  Chomiuk et al. 2011)—therefore we expect no radio emission if there is no more mass to interact with outside of R_w.

However, the late-time light curve of SN 2015bn could indicate continued interaction with CSM further out from the progenitor, for example a tail of lower density material causing the change in decline rate after around 70 days. In Figure 21, we showed a model that was powered by continued interaction with an extended outer region of CSM that had $\dot{M}/{v}_{10}=0.05$ M${}_{\odot }$ yr⁻¹, where v₁₀ is the wind velocity normalized to 10 km s⁻¹. Several other SLSNe also show a more gradual decline at late times (Figure 18). Radio observations at this phase can help to distinguish whether the slow luminosity decline is indeed powered by continued interaction, rather than ⁵⁶Ni-decay or a magnetar. If interaction is the culprit, as would be confirmed by bright radio emission, we can then use the radio luminosity to probe the mass-loss history of the progenitor.

We tested this scenario by comparing our deep radio limit to the predictions of interaction-powered models. Based on our estimate from Figure 21 (see also Table 3), we assume that the late-time luminosity is powered by continued interaction with a wind density profile located outside of ${R}_{{\rm{int}}}={10}^{15}$ cm. We can safely ignore the properties of any inner, denser regions of CSM, as the radio observations were carried out 200 days after the transition to the shallower decline rate, indicating the shock should long have left the dense region and will interact only with the outer material. Our model, which includes free–free absorption and synchrotron self-absorption (for further details on the model, see Kamble et al. 2014, 2016), conservatively assumes a shock velocity of 5000 km s⁻¹, which is at the lower end of the range we find in our fits to the optical light curve.

The radio SEDs of these models are shown in Figure 26. Comparing to our limits for SN 2015bn, we can see that it is our high-frequency observation that proves to be constraining. Therefore we compare the 22 GHz observation to the radio light curves in Figure 27. These models are calculated at 24 GHz to give a good match to the rest-frame frequency of the data. As demonstrated by the figure, the observational limit excludes models with mass-loss rates of ${10}^{-2.7}\lesssim \dot{M}/{v}_{10}\lesssim {10}^{-2}$ M${}_{\odot }$ yr⁻¹. Using a faster shock velocity does not greatly impact upon these results, but would generate brighter radio light curves that would rule out a slightly wider range in mass-loss rates. The largest mass-loss rate that we can exclude, $\dot{M}/{v}_{10}=0.01$ M${}_{\odot }$ yr⁻¹, is close to our estimate for that needed to power the late-time optical light curve through continued ejecta-CSM interaction. Models with $\dot{M}/{v}_{10}\approx 0.1$ M${}_{\odot }$ yr⁻¹ reach radio maximum at 10³ days after explosion. Therefore radio observations at $\nu \gtrsim 20\;{\rm{GHz}}$, carried out over the coming months and years, corresponding to progressively denser CSM, should be able to confirm or rule out the continued-interaction model we used to fit the optical light curve. This is perhaps the first easily testable prediction that can potentially distinguish between interaction- and magnetar-powered models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Inspection of Figure 26 shows that for high mass-loss rates ($\dot{M}/{v}_{10}\sim 0.1$ M${}_{\odot }$ yr⁻¹), the synchrotron flux is expected to peak in the millimeter range, suggesting that future observations of SLSNe should use both radio and millimeter observations to constrain the SED. The Atacama Large Millimeter/sub-millimeter Array (ALMA) should be able to reach the required sensitivity for the next SLSN at $z\approx 0.1$.