SPECTROPOLARIMETRY OF SUPERLUMINOUS SUPERNOVAE: INSIGHT INTO THEIR GEOMETRY

The last six years have seen the surprising discovery of new classes of intrinsically bright supernovae. They show absolute peak magnitudes of M ∼ −21 (e.g., Chomiuk et al. 2011; Quimby et al. 2011; Inserra & Smartt 2014) and a tendency to occur in dwarf, metal-poor galaxies (e.g., Chen et al. 2013, 2016; Lunnan et al. 2014; Leloudas et al. 2015b). They do not exhibit the typical narrow spectroscopic features of strongly interacting supernovae and they are usually referred to as superluminous supernovae (SLSNe; see the review of Gal-Yam 2012).

They can be divided into two main groups, hydrogen free—and hence labeled SLSNe I—and hydrogen rich—therefore called SLSNe II. The first class is better studied: they show blue continua at maximum light, a distinctive "W" feature due to O ii at early epochs and, at about 30 days after peak, they are spectroscopically similar to normal or broad-lined SNe Ic at peak luminosity (Pastorello et al. 2010); consequently, they are also usually labeled SLSNe Ic (e.g., Inserra et al. 2013). Additionally, SLSNe Ic show different light curve behavior and, as a consequence, have been divided in the subgroup of fast evolving (e.g., SN 2005ap, SN 2010gx; Pastorello et al. 2010; Quimby et al. 2011) and that of slow evolving (e.g., SN 2007bi, PTF12dam; Gal-Yam et al. 2009; Nicholl et al. 2013). The Type II SLSNe are fewer in number and they also show blue continua at maximum light, together with broad hydrogen features. At about 20 days after peak they show some resemblances to normal Type IIL, although by definition they are several magnitudes brighter (see Inserra et al. 2016, for a review).

Despite the increasing number of objects found every year, the nature of their explosion and the progenitor scenario are still debated. The favored scenario for all the types is that of an explosion driven by a magnetar as a central engine, which deposits energy into the supernova ejecta and significantly enhances the luminosity (Kasen & Bildsten 2010; Woosley 2010; Dessart et al. 2012). However, alternative scenarios such as the accretion onto a central black hole (Dexter & Kasen 2013), the interaction with a dense circumstellar medium (CSM; Chatzopoulos et al. 2012), and a pair-instability explosion are still feasible alternatives. The pair-instability mechanism is only physically plausible for some of the objects, in particular those with the slowest evolving light curve (e.g., see Gal-Yam et al. 2009; Kozyreva & Blinnikov 2015)

A powerful diagnostic for distinguishing between scenarios is polarimetry since it can unveil information on the geometry of the explosion and hence increase our understanding of the transients. Imaging polarimetry of a fast-evolving SLSNe Ic has been reported by Leloudas et al. (2015a), but no evidence of asymmetries was found. However, spectropolarimetry offers a more in-depth analysis of the geometry of SN explosions (see Wang & Wheeler 2008, for a review). Quantitative modeling of the lightcurves of SLSNe Ic indicates that the data are well matched by the explosion of massive progenitors with a central engine (e.g., as in Inserra et al. 2013; Nicholl et al. 2013). If this were true, we would expect to see an intrinsically asymmetric explosion characterized by a dominant polarization angle as observed for stripped-envelope SNe (e.g., Wang et al. 2001; Maund et al. 2007; Tanaka et al. 2012). Indeed many core-collapse SNe show large-scale departures from axisymmetry (e.g., Leonard et al. 2006; Wang & Wheeler 2008).

If a magnetar is the central engine that powers the extreme luminosity, then the strong magnetic field (B ∼ 10¹⁴ G) and rapid rotation could lead to asymmetries (Chen et al. 2016). Such asymmetries are potentially stronger than in normal stripped-envelope SNe. Detection of such signatures may suggest a core-collapse origin combined with an asymmetric magnetar energy injection process. On the other hand, alternative scenarios such as the interaction with a CSM would exhibit spectropolarimetric evolution similar to those of SNe IIn (e.g., Leonard et al. 2000; Wang et al. 2001). These signatures are different than that of stripped-envelope SNe, since they typically show unpolarized broad lines and loops across spectral lines arising from the CSM (Wang & Wheeler 2008).

An axisymmetric ejecta could be the consequence of an aspherical production of energy and momentum from an explosion due to magnetohydrodynamic jets (Khokhlov et al. 1999) or magnetoturbulence (e.g., Mösta et al. 2015), accretion flow around the central neutron star (e.g., Chevalier 1989), asymmetric neutrino emission (e.g., Wheeler & Akiyama 2010; Müller 2015), or a combination of these, or the fact that the material could be ejected in clumps.

In this paper we present our spectropolarimetric observations, as well as the interpretation of the geometry of the brightest known SLSN, namely SN 2015bn, which belongs to the group of slow-evolving SLSNe Ic.

SN 2015bn was discovered by the Catalina Real-time Transient Survey (CRTS, Drake et al. 2009 with IDs CSS141223-113342+004332 and MLS150211-113342+004333) on 2014 December 23 and subsequently by the Pan-STARRS Survey for Transients (PSST; with ID PS15ae, Huber et al. 2015) on 2015 February 15 UT. It was classified on the 2015 February 17 by the the Public ESO Spectroscopic Survey for Transient Objects (PESSTO; Smartt et al. 2015) as a Type Ic SLSN at z = 0.11 similar to several SLSNe Ic around maximum light (Le Guillou et al. 2015). The main PESSTO follow-up campaign is presented by Nicholl et al. (2016), where the redshift z = 0.1136 has been confirmed from the detection of narrow lines from the host galaxy. Adopting H₀ = 72 km s⁻¹, Ω_M = 0.27, and Ω_Λ = 0.73 as standard cosmological parameters, this corresponds to d_L = 514.3 Mpc. Alexander et al. (2016) reported that no radio continuum emission from the SN was detected by the Very Large Array on 2015 December 11.67 UT, with a 3σ upper limit of 40 μJy at a mean frequency of 21.85 GHz. This corresponds to a luminosity limit of of L_ν ≲ 1.3 × 10²⁸ erg s⁻¹ Hz⁻¹ at z = 0.1136, which sets a limit on the mass-loss rate ($\dot{M}\lt {10}^{-1.8}$ M_⊙ yr⁻¹) of the CSM in the case of the interaction scenario (Nicholl et al. 2016).

Spectropolarimetric observations of SN 2015bn were conducted with the Very Large Telescope (VLT) + the FOcal Reducer/low dispersion Spectrograph 2 (FORS2)¹ on the 2015 February 22 and April 20, which correspond to −23.7 and +27.5 days with respect to the maximum light in the g-band (MJD 57102.5, Nicholl et al. 2016). Both epochs are close to the broad peak of the SLSN light curve. A log of the observations is given in Table 1. We used a half-wave retarder plate at four angles (0°, 225, 45° and 675) through three iterations per epoch. All spectra were observed with the grism 300V and blocking filter GG453 in order to avoid second order contamination, giving an observed spectral range of 4400–9200 Å, and a dispersion of 1.68 Å pix⁻¹. A slit width of 1.0 arcsec resulted in a resolution of 13 Å. Data were reduced in a standard fashion (including trimming, overscan, bias correction, and flat-fielding) using standard routines within IRAF.²

Optimal extraction of the spectra was adopted. In addition, ordinary and extraordinary beams were processed separately. Wavelength calibration was performed using spectra of comparison lamps acquired with the same configurations as the SN observations. Atmospheric extinction correction was based on tabulated extinction coefficients for the telescope site. Flux calibration was performed using spectro-photometric standard stars observed on the same nights with the same setup as the SNe. The flux calibration was also checked by comparison with the photometry provided in Nicholl et al. (2016), integrating the spectral flux transmitted by standard griz filters—using the python program sms (Inserra et al. 2016)—and adjusted by a multiplicative factor when necessary. The resulting flux calibration matches the measured photometry to within 0.1 mag. Instrumental polarization and the position angle offset were checked by observing polarized and unpolarized standard stars during the second observation.

Stokes parameters were computed (Q and U) and error estimates were made using the procedures described by Patat & Romaniello (2006). The Stokes parameters were re-binned to 15 Å to increase the signal-to-noise ratio (S/N).

The Galactic reddening toward the position of the SN is E(B-V) = 0.02 (Schlafly & Finkbeiner 2011). No Na id λλ5890,5896 from the host galaxy was detected in our spectra, which is consistent with SLSNe (e.g., Inserra & Smartt 2014; Vreeswijk et al. 2014; Chen et al. 2016; Inserra et al. 2016) since they usually do not have reddening from the dwarf host galaxies. We will assume a total reddening equal only to the Galactic component. The assumption of both low host galaxy and Galactic reddening is supported by (i) the color of the spectra—T_BB ∼ 10,000 and ∼9400 K for the first and second epoch, respectively—which are similar to those of comparable SLSNe Ic at around the broad peak phase (Inserra et al. 2013; Nicholl et al. 2016), and (ii) the measurement of the host galaxy narrow emission lines of Hα, Hβ, and the resulting Balmer decrement, which points to E_Host(B – V) = 0 (Nicholl et al. 2016).

For the adopted reddening, the Serkowski–Whittet law allows us to place a limit of P_ISP < 0.2% (Serkowski et al. 1975; Whittet et al. 1992) on the interstellar polarization (ISP) due to intervening dust in the line of sight to SN 2015bn. In addition, we found three stars from the Heiles (2000) catalog within three degrees of the SN position, which have P = 0.07%, P = 0.10%, and P = 0.05%, corresponding to a mean of $\bar{{\text{}}P}=0.07\pm 0.04 \%$, in full agreement with the low Galactic ISP along the line of view.

Furthermore, as an independent test, we have derived the ISP value directly from our spectropolarimetric observations, assuming that the intrinsic polarization in regions of strong emission lines or line blanketing is negligible (Trammell et al. 1993). With this assumption, polarization signals in these spectral regions can be attributed to the ISP. We identified the region of the emission of Si ii λ6355 and the region around 4500 Å affected by Fe ii and Fe iii blending (see Section 5) as suitable for such analysis. The Stokes parameters were calculated over a range of 100 Å in the binned (15 Å) regions of 4670–4770 Å and 6300–6400 Å. In the first epoch we measured Q_ISP = −0.15% and U_ISP = −0.06%, which gives us a polarization level P_ISP = 0.16% consistent with the previous limits. For our analysis, we applied this correction to the spectra of both the first and the second epoch—as the latter does not present the iron blanketing.

To evaluate the errors on our polarimetry, we measured the rms of the normalized flux differences of our ordinary and extraordinary beams (see Equations (4) and  (5) in Patat & Romaniello 2006) and then propagate those to the Stokes parameters and the polarization vector. Since in our data set η = P₀(S/N) ∼ 8, where P₀ is the input polarization and S/N is the signal-to-noise ratio (Sparks & Axon 1999; Patat & Romaniello 2006), we can check the validity of our error estimates by comparing them to the analytic expression for the absolute error in P given by (Patat & Romaniello 2006, see their Equation (10)). We find that the analytic error estimates are in good agreement with our error calculations, which we adopt throughout the following. Furthermore, due to the large η value, we can safely assume that both P and σ_P follow a Gaussian distribution and do not require debiasing corrections (Sparks & Axon 1999).

The spectroscopic evolution of SLSNe Ic has been well sampled (e.g., Pastorello et al. 2010; Inserra et al. 2013), especially that of slow-fading ones (e.g., Nicholl et al. 2013, 2016; Vreeswijk et al. 2014) with the evolution from pre- to post-peak spectra mainly showing the disappearance of the characteristic O ii lines soon after peak and the decrease of ionization state of the metal lines. As already shown in Nicholl et al. (2016), SN 2015bn follows such evolution.

Our first spectrum was taken 24 days before peak magnitude in g-band. A blackbody fit to the available spectral range yields T ∼ 10,000 K, which corresponds to a radius of 6.5 × 10¹⁵ cm and blackbody peak bluer (2897.8 Å) than covered by our wavelength range. The spectrum is dominated by metal lines: the Ca H&K lines are visible on the blue edge of our spectrum, redward of which are two double peak absorption features at 4400 and 5000 Å, most likely due to a combination of Fe iii and Fe ii with a possible contribution of O ii, which could be excited by non-thermal processes (Mazzali et al. 2016). Redward, the spectrum is again dominated by iron up to 5500 Å, while Na id and Si ii λ6355 are visible in the center part around 6000 Å. The region around 7500 Å shows a shallow O i absorption line, together with another shallow P-Cygni profile at 7200 Å present only in the slow-fading SLSNe Ic and tentatively identified as C i (see Figure 11 in Nicholl et al. 2016). Following the modeling of SLSNe Ic spectra reported by Mazzali et al. (2016), an alternative identification for the 7200 Å line could be O ii.

The second spectrum—taken 28 days after the maximum light—has a slightly lower blackbody temperature (T ∼ 9400 K) and larger blackbody radius 8 × 10¹⁵ cm compared to the first. At this epoch Fe iii has been replaced by Fe ii (see top panels of Figure 1). The noticeable changes are in the region blueward of Si ii. Mg i] λ4571 appeared in emission and replaced Fe iii and—together with the increase of Fe ii features with respect to the previous epoch—has reshaped the profile of that region that no longer shows the two absorption features. The other distinguishable change is the shift to the red of the broad profile around 5300 Å that is due again to Fe ii lines replacing their higher ionized ions. The other two P-Cygni profiles redder than Si ii are still present and due to the same elements of the first epochs, even if the feature around 7200 Å could now have a small contribution from [Ca ii] λλ7291, 7323 as also shown by Nicholl et al. (2016).

Zoom In Zoom Out Reset image size

Our spectropolarimetric observations of SN 2015bn, not corrected for the ISP, are shown in Figure 1. In both our spectra we detect clear polarization signals. In the first epoch, the typical degree of polarization increases from ∼0.5% in the blue to a maximum of P ∼ 1.2% at ∼7000–7500 Å. However, in the post-maximum spectrum we find the opposite behavior: a decrease of the polarization from the blue (P ∼ 1.4%) to the red (P ∼ 0.5%).

Most of the peaks in the polarization signal are roughly coincident with the absorption minima of identifiable P-Cygni profiles in the spectra: at the first epoch, there is a clear increase in polarization across the Fe iii absorption feature at about 5100 Å and also across both the oxygen lines at 4400 Å and 7700 Å . At the second, post-maximum epoch (+27.5 day), the degree of polarization across O i (P ∼ 0.9%) remains similar to that found at the first epoch (P ∼ 1.2%). However, there are significant changes to all the features blueward of O i. The Si ii line, which shows a more prominent P-Cygni profile in the +27.5 day flux spectrum than the −23.7 day one, reaches a peak polarization signal of ∼1%, while the Fe ii 5000 Å features reaches P ∼ 1.7%.

The high-polarization feature at about 7000 Å corresponds to the blue end of a shallow absorption trough in both the spectra, which was tentatively identified as carbon by Nicholl et al. (2016) or it could alternatively be associated with oxygen (see Section 5).

The general evolution from pre-maximum to post-maximum phase shows an increase of polarization—especially for the metal lines at blue wavelengths—as usually observed in stripped-envelope SNe (e.g., Wang et al. 2003; Maund et al. 2007; Reilly et al. 2016).

6.1. Dominant Axis and Asymmetry

In Figure 2 we show our spectropolarimetric data set in the Q – U plane. As before, we note that the polarization reaches a maximum of about 1.4% in the first epoch and almost 2% in the second.

Zoom In Zoom Out Reset image size

The distribution of our data points in the Q – U plane clearly suggests a preferred direction, which we identify by fitting the data with a least-squares fit weighted by the observational errors following the procedure of Wang et al. (2003). The best fits to the available data set give a dominant axis with a position angle (θ_d = 1/2 × arctanβ) of −8° ± 1 for the pre-maximum and −6° ± 1 for the post-maximum epoch with respect to a north–south axis on the night sky.³ We note that, in the first epoch, the dominant axis almost crosses the origin and passes through the ISP measurement (black star in Figure 2), suggesting that our estimate for the low level of ISP was correct. On the other hand, the best fit of the second epoch is slightly different than the previous and does not cross the ISP. In principle, a difference might indicate an intrinsic change of the dominant axis from the first to the second epoch and hence a different dominant orientation between the external and inner layers of the ejecta. Alternatively, the difference between the fits can be due to the larger uncertainties on Q and U in the second epoch, especially in the red part of the spectrum. However, the change in the angle of the best fit between the two epochs is sufficiently small, and comparable with the errors, so it can be attributed to measurement uncertainties only. We note that if we force our fit to pass through the ISP in the second epoch, we obtain an angle of −11° ± 1, which is still comparable to that of the first epoch.

We note that our analysis clearly suggests that the data at the first epoch cross the origin in the Q – U plane (from second to fourth quadrant). As discussed by Patat et al. (2010), this often suggests that the ISP may have been inaccurately estimated and that the true value of the ISP should be at one of the two extremes of the data distribution in the Q – U plane. In our case, however, that would imply a 1% < P_ISP < 2%, which is inconsistent with the low reddening observed toward the SN (see Section 4). Thus we favor an interpretation in which the data distribution in the Q – U plane, including the crossing between quadrants, is intrinsic to the ejecta and not a consequence of significantly underestimated ISP.

It is convenient to project the polarization vectors in a new coordinate system that has a component along the dominant axis

$$\begin{eqnarray}&&{Q}_{\mathrm{rot}}=(Q-{Q}_{\mathrm{ISP}})\cos (2\theta )+(U-{U}_{\mathrm{ISP}})\sin (2\theta ),\end{eqnarray} \tag{ 1 }$$

which represents global geometric deviations from spherical symmetry, and an orthogonal component to the dominant axis

$$\begin{eqnarray}&&{U}_{\mathrm{rot}}=-(Q-{Q}_{\mathrm{ISP}})\sin (2\theta )+(U-{U}_{\mathrm{ISP}})\cos (2\theta ),\end{eqnarray} \tag{ 2 }$$

which shows physical deviations from the dominant axis (for further details see Leonard et al. 2001; Wang & Wheeler 2008). In Figure 3 we show our polarization vectors in the new system. For both epochs we used a single angle of θ = 8°. This makes clear that, while our data show a strong deviation from spherical symmetry (seen in Q_rot), there is only weak evidence of physical deviation around the dominant axis (i.e., U_rot is fairly consistent with zero, suggesting that the ejecta are approximately axisymmetric). At both epochs, Q_rot shows a wavelength dependence decreasing from blue to red. At the first epoch, we find a clear sign reversal with Q_rot changing from positive values for blue wavelengths (≲6000 Å) to negative values in red. The second epoch is more polarized than the first and does not show a sign reversal.

Zoom In Zoom Out Reset image size

SN 2015bn would be classified as spectropolarimetry type D0 in the terminology of Wang & Wheeler (2008) since it shows a distinct dominant axis and a distribution orthogonal to that consistent with observational noise. In addition, the stronger lines do not exhibit a loop in the Q – U diagram, hence they do not show large significant changes in amplitude and position angle.

6.2. Comparison with LSQ14mo

6.3. Modeling

As discussed in the previous section, polarization spectra for SN 2015bn are characterized by a well-defined axis of symmetry. In addition, the polarization component along the dominant axis, Q_rot, shows a strong wavelength dependence at both epochs and a pronounced evolution between pre- and post-maximum. In the following we present a simple toy model that can account for such observed polarization signatures. We note that this model is only intended to characterize the observed data and that the results presented here may not be unique: similar results may be obtained with different geometries and compositions. Our primary goal is to constrain the physical mechanisms responsible for (i) the wavelength dependence and (ii) the time evolution of the overall pseudo-continuum polarization level, rather than to model polarization features associated with individual spectral transitions.

To interpret the polarization spectra for SN 2015bn, we will consider simple ellipsoidal geometries using the Monte Carlo toy-code presented by Bulla et al. (2015). Monte Carlo packets, representing bundles of photons, are created unpolarized at a spherical photosphere of radius R_ph, emitted assuming constant surface brightness and propagated throughout a prolate ellipsoidal envelope defined by

$$\begin{eqnarray}&&\displaystyle \frac{{x}^{2}}{{a}^{2}}+\displaystyle \frac{{y}^{2}}{{b}^{2}}+\displaystyle \frac{{z}^{2}}{{c}^{2}}=1,\,a=b\lt c.\end{eqnarray} \tag{ 3 }$$

The latter can be rewritten in cylindrical coordinates as

$$\begin{eqnarray}&&\displaystyle \frac{{r}^{2}}{{A}^{2}}+{z}^{2}={c}^{2},\end{eqnarray} \tag{ 4 }$$

where we have introduced the axis ratio A = a/c < 1. Our calculations adopt ellipsoidal isodensity surfaces

$$\begin{eqnarray}&&\rho (\xi ,t)=\displaystyle \frac{{\rho }_{0}(t)}{{\xi }^{n}},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&\xi =\sqrt{\displaystyle \frac{{r}^{2}}{{A}^{2}}+{z}^{2}}.\end{eqnarray} \tag{ 6 }$$

For all our calculations we use a power-law index n = 7 as suggested by previous studies (e.g., Mazzali et al. 2016). On their journey throughout the ejecta, Monte Carlo packets can interact with matter and change their polarization state. Two sources of opacity are used: a gray electron scattering opacity and a wavelength-dependent line opacity (see below). The former is assumed to partially polarize the radiation, while the latter is regarded as a depolarizing contribution (see e.g., Jeffery 1989, 1991). Polarization spectra are extracted for an equatorial viewing angle (that is along the semi-minor axis of the prolate ellipsoid). Moving the observer orientation from the equatorial plane toward the pole would result in smaller polarization levels and thus would require larger asphericities than those predicted below.

The electron scattering coefficient and line opacity distribution adopted in our calculations are estimated using the one-dimensional tardis radiative transfer code (Kerzendorf & Sim 2014). In particular, calculations with two different ejecta compositions have been carried out: one (hereafter referred to as Model A) adopting the ejecta composition proposed by Mazzali et al. (2016) and one (Model B) assuming an extreme case of ⁵⁶Ni-dominated ejecta at the time of explosion. In the latter case, we derive the relative abundances of each element by calculating the fraction of Ni that has decayed to Fe and Co at both epochs of SN 2015bn. The specific ejecta compositions used for each calculation are reported in Table 2. tardis provides values for the electron scattering coefficient and the Sobolev optical depth of each transition in each zone from which tables of wavelength-dependent expansion opacities can be calculated (Karp et al. 1977; Friend & Castor 1983; Eastman & Pinto 1993). The opacities obtained for tardis zones close to the photosphere in tardis are plotted in Figure 4. At both epochs, Model A produces only a weak wavelength dependence in the line opacity distribution across the wavelength interval covered by our observations (∼4000–8000 Å), with values that are a factor ∼10 lower than the electron scattering opacity (${k}_{\mathrm{sc}}^{A,1}=0.047\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ and ${k}_{\mathrm{sc}}^{A,2}=0.027\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ at the first and second epoch, respectively). In contrast, the line opacity distribution extracted for Model B is strongly wavelength dependent in this region. While spectral regions around 8000 Å are still dominated by electron scattering opacity (${k}_{\mathrm{sc}}^{B,1}=0.015\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ and ${k}_{\mathrm{sc}}^{B,2}=0.014\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$), the line opacity contribution steadily increases toward the bluer regions of the spectrum and becomes dominant around 4000 Å.

Zoom In Zoom Out Reset image size

Table 2.  Parameters of the Models for SN2015bn

Model	Epoch	Log10 Lbol	vph	vout	ρph	He	C	O	Ne	Mg	Si	S	Ca	Ti	Fe	Co	Ni
		(erg s−1)	(km s−1)	(km s−1)	(g cm−3)
A	02/22	44.2	8e8	1.6e9	7.3e-14	0.10	0.40	0.475	0.02	7e-4	0.002	5e-4	5e-5	2e-5	5e-4	3e-4	1e-5
A	04/20	44.1	8e8	1.6e9	3.6e-14	0.10	0.40	0.475	0.02	5e-4	0.002	5e-4	5e-5	2e-5	5e-4	3e-4	1e-5
B	02/22	44.2	8e8	1.6e9	7.3e-14	0	0	0	0	0	0	0	0	0	0.414	0.585	3.5e-4
B	04/20	44.1	8e8	1.6e9	3.6e-14	0	0	0	0	0	0	0	0	0	0.626	0.374	1e-6

Download table as:  ASCIITypeset image

Figure 5 reports results of two simulations relative to the first epoch of SN 2015bn. Both calculations adopt opacities from Model A and assume an ellipsoidal envelope with axis ratio A = 0.88 and outer boundary at ξ_out = 2 × R_ph. This value of the axis ratio was chosen to match the continuum polarization around 7000 Å (see below). The maximum electron scattering opacity to the boundary

$$\begin{eqnarray}&&{\tau }_{\max }={\int }_{{R}_{\mathrm{ph}}}^{{\xi }_{\mathrm{out}}}{k}_{\mathrm{sc}}\,\rho (\xi ,t){dz}\end{eqnarray} \tag{ 7 }$$

is set to unity. Ejecta are modeled 70 days after explosion and the photosphere placed at ${v}_{\mathrm{ph}}^{1\mathrm{st}}=8\times {10}^{3}$ km s⁻¹ (Nicholl et al. 2016). The density at ξ = R_ph and the mass in the envelope are ${\rho }^{1\mathrm{st}}=2.7\times {10}^{-14}$ g cm⁻³ and ${M}^{1\mathrm{st}}=2.5$ M_⊙, respectively. The first polarization spectrum (red dashed line) corresponds to a calculation including only electron scattering (${k}_{\mathrm{sc}}^{A,1}$) and provides a good match to the polarization level detected for SN 2015bn in the wavelength region around 7000 Å that is devoid of strong optical features. The second polarization spectrum (solid blue line) shows instead a simulation in which both electron scattering and line opacities are included. Specifically, we adopt a polynomial fit to the line opacity distribution as representative of the pseudo-continuum absorption component (see upper left panel of Figure 4). Compared to the line-free calculation, the latter produces lower polarization levels since the line opacity acts as a depolarizing contribution. This effect becomes stronger moving from redder to bluer wavelengths, as a consequence of the increasing contribution of line to the total opacity (see upper left panel of Figure 4). The match between predicted and observed polarization levels, however, is poor as the wavelength dependence in the line opacity of Model A is not sufficiently strong to bring the overall Q_rot level toward positive values below ∼5500 Å, as observed for SN 2015bn.

Zoom In Zoom Out Reset image size

Compared to Model A, line opacities extracted for Model B are characterized by a more pronounced wavelength dependence (see right panels of Figure 4). Therefore, we construct a new model that includes opacities from both Model A and Model B. Specifically, we divide the envelope into two different ellispoidal zones:

• 1.  
  A "metal-rich" inner region modeled as a prolate ellipsoid with ${\xi }_{\mathrm{in}}=1.5\times {R}_{\mathrm{ph}}$. Opacities in this region are selected from Model B. Specifically, we adopt ${k}_{\mathrm{sc}}^{B,1\mathrm{st}}$ (pre-maximum) and ${k}_{\mathrm{sc}}^{B,2\mathrm{nd}}$ (post-maximum) as electron scattering coefficients and use a polynomial fit to the line opacity distribution for the pseudo-continuum absorption component (see right panels of Figure 4).
• 2.  
  A "metal-poor" outer region modeled as a prolate ellipsoid with ${\xi }_{\mathrm{out}}=2\times {R}_{\mathrm{ph}}$. Opacities in this region are selected from Model A. Here, we adopt electron scattering coefficients ${k}_{\mathrm{sc}}^{A,1\mathrm{st}}$ (pre-maximum) and ${k}_{\mathrm{sc}}^{A,2\mathrm{nd}}$ (post-maximum) and a polynomial fit to the line opacity distribution for the pseudo-continuum absorption component (see the left panels of Figure 4).

The shapes of the two zones (i.e., the axis ratios A_in and A_out) are taken as free parameters in our models, together with the maximum electron scattering opacity to the boundary τ_max. A suitable choice for τ_max is selected (see below) and the axis ratios A_in and A_out chosen to reproduce the observed polarization spectra.

The left panel of Figure 6 shows the results of our simulations compared to the first epoch of SN 2015bn. As in the one-zone model presented above, these calculations are carried out at 70 days after explosion and assume ${v}_{\mathrm{ph}}^{1\mathrm{st}}=8\times {10}^{3}$ km s⁻¹ and ${\tau }_{\max }^{1\mathrm{st}}=1$. The inner and the outer ellipsoidal envelopes have the same axis ratio, namely ${A}_{\mathrm{out}}^{1\mathrm{st}}={A}_{\mathrm{in}}^{1\mathrm{st}}=0.88$. Estimated values for the density at ξ = R_ph and for the mass in the envelope are ${\rho }^{1\mathrm{st}}=7.3\times {10}^{-14}$ g cm⁻³ and ${M}^{1\mathrm{st}}=6.7$ M_⊙, respectively. The first simulation (dashed red line) includes only electron scattering in the ejecta and is found to reproduce the degree of polarization observed for SN 2015bn around 7000 Å (Q_rot ∼ −1%). While the gray electron scattering opacity assumed in the latter calculation yields a constant degree of polarization as a function of wavelength, including contributions from line absorption (solid blue line) results in a strong wavelength dependence of the pseudo-continuum level. This causes a sign reversal of ${Q}_{\mathrm{rot}}$ around 5000 Å, in good agreement with what is observed for SN 2015bn. The similarity between our predictions and the overall pseudo-continuum level of SN 2015bn thus suggests that an increasing contribution of line opacities toward the bluer regions of the spectra is responsible for the strong wavelength dependence observed.

Zoom In Zoom Out Reset image size

The right panel of Figure 6 shows the results of our simulations relative to the second epoch of SN 2015bn. The latter is modeled 121 days after explosion and with the photosphere placed at ${v}_{\mathrm{ph}}^{2\mathrm{nd}}=7\times {10}^{3}$ km s⁻¹ (Nicholl et al. 2016). Following the ejecta expansion, the density at ξ = R_ph and the maximum electron scattering opacity drop to ${\rho }^{2\mathrm{nd}}=3.6\times {10}^{-14}$ g cm⁻³ and ${\tau }_{\max }^{2\mathrm{nd}}=0.61$, respectively, while the mass in the envelope increases to ${M}^{2\mathrm{nd}}=11.9$ M_⊙. As shown in Figure 6, keeping the axis ratios of the "metal-rich" and "metal-poor" regions fixed to the values they have at the first epoch leads to polarization levels that are clearly inconsistent with those observed in SN 2015bn (dot-dashed green line). If we instead keep the same axis ratio for the outer region (${A}_{\mathrm{out}}^{2\mathrm{nd}}={A}_{\mathrm{out}}^{1\mathrm{st}}=0.88$) and assume a more aspherical inner region (${A}_{\mathrm{in}}^{2\mathrm{nd}}=0.60$), we obtain a polarization spectrum in good agreement with data (solid blue line).

To summarize, our ellipsoidal toy model suggests that (i) the strong wavelength dependence observed in both polarization spectra of SN 2015bn is given by a strong increase in the line opacity from redder to bluer regions of the spectra, which has been previously predicted for models of SNe II (Dessart & Hillier 2011), and (ii) the time evolution of the polarization level from pre- to post-maximum is given by a change in the asphericity of the inner layers (from ${A}_{\mathrm{in}}^{1\mathrm{st}}=0.88$ to ${A}_{\mathrm{in}}^{2\mathrm{nd}}=0.60$). The ejecta geometries proposed by our calculations are reported in Figure 7. The implications of these results for the different progenitor scenarios will be discussed in Section 7.

Zoom In Zoom Out Reset image size

We stress again that the simple ellipsoidal modeling developed here is only intended to broadly characterize the data and frame their interpretation. In particular, we did not attempt to model any discrete spectral features nor investigate the extent to which the two epochs are consistent with evolution in homologous expansion. Instead we have considered the two epochs individually.

Several scenarios have been proposed to account for the exceptional luminosity of SLSNe. Among these, the most successful for reproducing light curves, as well as velocity and temperature evolution, has been that of a central engine, likely a magnetar, that deposits its energy into the supernova explosion (Kasen & Bildsten 2010; Woosley 2010; Inserra et al. 2013; Nicholl et al. 2013). However, other scenarios including the interaction of the ejecta with CSM due to a pulsational pair-instability SN (PPISNe, e.g., Woosley et al. 2007) or a black hole as the central engine (Dexter & Kasen 2013) cannot be ruled out in explaining the variety of SLSNe observables. Our spectropolarimetric observations provide a new insight into the geometry of an SLSN explosion and we can now consider the implications of this for identifying the explosion mechanism responsible for their extreme luminosity.

The data show a deviation from spherical symmetry and the presence of a dominant axis, which implies an axisymmetric configuration of the ejecta. Furthermore, the U_rot polarization vector in the new rotated system—where the Q_rot vector is along the dominant axis and the U_rot is orthogonal to the previous one (see Section 6.1)—indicates that there is almost no physical deviation from the dominant axis. To explain the observed increase in the Q_rot vector from the pre-peak to the post-peak epoch (see Figure 3) we need an increase in the asphericity of the inner layers of the ejecta or the photosphere to recede into inner layers that are less spherical. The geometric parameters inferred from our toy ellipsoidal modeling of the inner zone (A${}_{\mathrm{in}}^{1}=0.88$ and ${A}_{\mathrm{in}}^{2}=0.60$; see Section 6.3) are certainly model dependent and so can only be taken as indicative: they depend in detail on the composition/line opacity adopted and on all parameters controlling the shape and density distribution. To fully understand degeneracies among these parameters, and to construct a self-consistent model that accounts for both epochs is beyond the scope of the simplistic approach used here. Nevertheless, the calculations illustrate that moderate departures from sphericity combined with substantial wavelength-dependent opacity can qualitatively account for the observed spectropolarimetric properties. Although we used ellipsoidal toy models to characterize the data (shown in Figure 7), it is important to note that alternative axisymmetric configurations, including bipolar ejecta morphologies, would also be able to match the observations.

Together, the data and modeling presented above suggest that, as time passes, the spectral forming region recedes into increasingly aspherical (but axisymmetric) ejecta layers. In the following, we consider the ramifications of this conclusions for our understanding of SLSN explosion models, bearing in mind that models involving a dominant axis are favored by the data.

Multidimensional simulations of magneto-rotational core-collapse SNe have already shown large-scale asymmetries (Mösta et al. 2014). In the case of a magnetar as the inner engine, the overall axisymmetry and the increase of asphericity as the photosphere recedes through the inner layers with time could be the consequence of several factors.

• 1.  
  It could be a consequence of a jet-like flow stalled within the core in a similar fashion to that proposed for SN 2008D/XRF 080109 (Maund et al. 2009). This is an event linked to X-ray flashes, which may also be related to γ-ray bursts but with a softer radiation outburst. Hypothetically, this suggestion could also be supported by the possible link between SLSNe Ic and ultra-long γ-ray burst (Greiner et al. 2015; Metzger et al. 2015; Kann et al. 2016).
• 2.  
  As the ejecta expands, its density decreases and the non-thermal radiation from the magnetar wind nebula can ionize more of the ejecta (Metzger et al. 2014). A magnetar wind nebula that is asymmetric could enhance such behavior.
• 3.  
  Two-dimensional simulations of magnetar-powered SLSNe (Chen et al. 2016) have shown that fluid instabilities and mixing occurs if energy is deposited in a relatively small mass of material, which could be, approximately, the mass of the inner ejecta. The amount of energy deposited needs to be comparable to that of the ejecta kinetic energy of the explosion. These instabilities and mixing could alter the geometry of the ejecta, giving an axysimmetric geometry reminiscent to that of the Crab Nebula, if the initial magnetar period is less than 3 ms (Chen et al. 2016). Notably, the magnetar models proposed for the bolometric light curve of SN 2015bn do satisfy this condition (2.1 and 1.7 ms for the two best fit shown in Nicholl et al. 2016).

In the case of a black hole as a central engine, the first scenario described above (the jet-like flow) would still be applicable. To date, there have been no multidimensional simulations of this mechanism and little can be added about the geometry of the ejecta. However, the inner engine mechanism would need to explain the wavelength dependence observed in SN 2015bn data.

Recently Kozyreva & Blinnikov (2015) have suggested that pair-instability SNe (PISNe) could still be a viable scenario to explain slow-evolving SLSNe Ic like SN 2015bn. Ejecta from such explosions have metal-rich inner ejecta. The attempt to model the wavelength dependence in our spectropolarimetry data suggests a strong line opacity, achievable with a very metal-rich composition. However, there is some tension between this ejecta composition and the observed spectra since they do not show very obvious signatures of strong line blanketing. This suggests that detailed modeling, including self-consistent spectra and polarization, is needed. Additionally, it is hard for such massive explosions to break the original spherical symmetry as shown by several multidimensional simulations (Fryer et al. 2001; Joggerst & Whalen 2011; Baranov et al. 2013; Chen et al. 2014a). A multidimensional simulation including rotation was actually able to produce axisymmetric inner ejecta for a highly rotating CO core (Chatzopoulos et al. 2013), but this was at the cost of lower luminosity and redder spectra than those of non-rotating PISNe. We note that one object showing similar behavior to that of highly rotating PISNe has been reported (PS1-14bj, Lunnan et al. 2016). However, a PISN explosion for SN 2015bn is hard to reconcile with all the data we have to date (see also Jerkstrand et al. 2016).

The last scenario to consider is that of the CSM interaction. As mentioned in Section 1, the most promising mechanism to explain the high luminosity in the framework of interaction is that of a pulsational pair-instability SN (PPISN), where thermonuclear outbursts due to a recurring pair-instability release massive shells of material that eventually collide with each other. The undulations observed in the SN 2015bn light curve (Nicholl et al. 2016), as in other slow-evolving SLSNe Ic (C. Inserra et al. 2016, in preparation; Lunnan et al. 2016), could be a consequence of collisions between multiple shells. Such collisions could lead to a strong mixing and Rayleigh–Taylor instabilities, creating a structure that  deviates from spherical symmetry behind the shock generated in the last collision (Chen et al. 2014b). However, such instabilities would not preserve an axis. Additionally, the outermost collision is like that between SN ejecta and a dense shell expelled by the star prior to its death as in the case of Type IIn. This was suggested for the case of SLSN IIn SN 2006gy (Woosley et al. 2007), which shows narrow H lines in the spectra throughout its evolution. Hence, it is difficult to reconcile the flux and polarization spectra with those of interacting SNe. Indeed, even if we assume that in our H- and He-free environment we should not observe narrow lines belonging to the unshocked external material. We should have observed an overall polarized continuum as in the other interacting supernovae (e.g., Patat et al. 2011) and a line polarization similar to those observed in normal, non-interacting SNe, which show an axisymmetric photosphere below the interaction. As reported in Sections 5 and 6 and shown in Figure 1 we do not observe such behavior.

We have presented the first spectropolarimetric data for an SLSN. We have gathered two epochs, one pre-peak at −23.7 days and the second 27.5 days after maximum light in the rest-frame.

Our analysis of these data indicates:

• 1.  
  the presence of a dominant axis and no physical departure from it;
• 2.  
  a strong wavelength dependence of the polarization spectra; and
• 3.  
  an increase in the mean degree of polarization from the first to the second epoch of observations.

We used our Monte Carlo toy-code to compare with the data in order to interpret the wavelength dependence and the time evolution of the polarization level. We calculated polarization spectra for prolate ellipsoidal geometries accounting for electron scattering and resonant line scattering. We were able to reproduce the wavelength dependence of the polarization by adopting a line opacity distribution for the inner ejecta that is rich in iron-group elements. This poses some challenges for the observed flux spectra, but further comprehensive modeling with full radiative transfer is required. We have also show that the evolution of the overall pseudo-continuum level can be replicated via two-zone aspherical ejecta models with an outer zone having A_out = 0.88 at both epochs and an inner zone showing an increase in the asymmetry from the first (A_in = 0.88) to the second (A_in = 0.60) epoch. These values are calculated for an equatorial viewing angle and thus should be considered as lower limits on the asphericities of the ejecta. Orientations away from the equatorial plane would lead to lower polarization signals and would then require smaller axis ratio values. In addition, we caution that the specific values derived here (A_in = 0.88–0.60 and A_out = 0.88) are model dependent. However, our findings qualitatively demonstrate that the SN 2015bn spectropolarimetry can be reproduced by an ellipsoidal geometry—or alternatively a bipolar geometry—with an inner region that increases its asymmetry as time passes.

The general trend for an increase in polarization with time (and the implied geometry) is reminiscent of those observed in stripped-envelope core-collapse SNe connected with γ-ray bursts, supporting the possible link between SLSNe Ic and ultra-long γ-ray bursts. Among all the suggested scenarios for SLSNe I, these new data tend to favor a core-collapse explosion and a central inner engine as an explanation for the inferred axisymmetric geometry, as well as the increase of the inner asymmetry. This could be achieved with an explosion of a rotating stellar progenitor. Currently, this scenario does not comfortably reproduce the wavelength dependence we observe, since it seems challenging to obtain enough line opacity. On the other hand, we did not explore other opacity options than a simple iron core and hence a more complete modeling is needed (including self-consistent spectra and polarimetry).

Despite the uncertainties that still exist with our relatively simple modeling approach, the data should drive new initiatives in how to model these data with hydrodynamic simulations. These models will address questions such as the the geometry (bipolar or ellipsoidal shapes) and how an engine-driven or interaction scenario could explain the wavelength dependence we observe (which we interpret as due to line opacity). We also need to know if other SLSNe I and broad-line SLSNe II show similar properties. The next step is to collect further observations of the future nearby (z < 0.2) SLSNe at least at two epochs, ideally one before and one post-maximum light. Predicting polarization signatures for multidimensional hydrodynamic explosion models will then be required to better investigate the geometry of these explosions.

This paper is based on observations made with ESO Telescopes at the Paranal Observatory under program ID 095.D-0611(A). C.I. thanks the ESO staff of the VLT, and especially Joe Anderson, for his competent support of this project in service mode and to have carried out these great observations. C.I. also thanks Matt Nicholl for kindly having provided information about SN 2015bn before the submission of his paper and for comments on SN 2015bn. C.I. and M.B. thank Emma Reilly, Ferdinando Patat, and Stefan Taubenberger for their useful discussions on the reduction and errors evaluation of spectropolarimetric data. C.I. also thanks Ting-Wan Chen and Anders Jerkstrand for useful insights on the topic of SLSNe. We thank PESSTO (the Public ESO Spectroscopic Survey for Transient Objects Survey), ESO programs 188.D-3003, 191.D-0935 for the prompt classification of the object. S.J.S. acknowledges funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement No. [291222] and STFC grants ST/I001123/1 and ST/L000709/1.