© ESO 2020

1. Introduction

Type IIn supernovae (SNe IIn) were introduced as a class by Schlegel (1990), who noted the presence of narrow hydrogen emission lines superposed on blue spectra. The emission lines that are characteristic of SNe IIn are generated from the interaction between rapidly expanding SN ejecta and a hydrogen-rich circumstellar medium (CSM), which drives the formation of a forward shock into the CSM and a reverse shock into the ejecta (see, e.g., Chevalier & Fransson 1994). Multiple components can form the emission lines. These include a narrow (v_FWHM ∼ 100 km s⁻¹) component originating from unshocked CSM, ionized by the high-energy photons emitted from the shock region. In addition, an intermediate component (v_FWHM ∼ 1000 km s⁻¹) has been observed to arise from the region between the forward shock and the reverse shock. A broad component (v_FWHM ∼ 10 000 km s⁻¹) can arise from the ionized expanding SN ejecta, or it is more likely due to electron-scattering broadening of the narrower component by an optically thick CSM (e.g., Chugai 2001; Huang & Chevalier 2018).

The picture is further complicated by the distribution of the CSM and the presense of dust. For example, the CSM could be aspherical and/or asymmetric (e.g., Stritzinger et al. 2012), while the presence of dust (e.g., Fox et al. 2011) in the ejecta and/or in the CSM could modify the emission line profiles. Dust also produces an excess of red flux as it redistributes bluer flux to near-infrared (NIR) and mid-infrared (MIR) wavelengths (see, e.g., Gall et al. 2014).

A variety of light-curve shapes and luminosity levels characterize the various SN IIn populations (see review in e.g., Smith 2017). In Taddia et al. (2013, 2015), three main subgroups are schematized as namely long-lasting SNe IIn or 1988Z-like objects, fast-declining SNe IIn-L or 1998S-like objects, and plateau SNe IIn-P or 1994W-like objects. Nyholm et al. (2020) identify two groups of slow and fast rising SNe IIn in their Palomar Transient Factory sample. The large range of luminosities goes from faint (M_r ∼ −16 mag) SNe IIn to superluminous type II SNe (SLSNe II) that appear to be extremely bright SNe IIn with M_r ∼ −22 mag (Nyholm et al. 2020). SN IIn light curves can be used, together with information from the spectral lines, such as CSM velocity, shock velocity, and ejecta velocity, to estimate the mass-loss rate, the ejecta mass, and the explosion energy of the SN progenitor that produces the CSM via its stellar winds and/or episodic mass-loss events (e.g., Moriya et al. 2013, 2014).

Taddia et al. (2015) find that different SN IIn subgroups seem to occur in different environments, with long-lasting SNe IIn occurring at lower metallicity than fast-declining SNe IIn. Long-lasting SNe IIn have similar environmental properties as those of SN impostors that are believed to be produced by luminous-blue-variable (LBV) eruptions. Transients likely arising from significant stellar eruptions sometimes occur prior to the birth of SNe IIn. Examples of such objects include the outbursts observed prior to and at the position of SN 2006jc (Pastorello et al. 2007) and several SNe IIn studied by PTF/iPTF (Ofek et al. 2014a; Nyholm et al. 2017). SN 2009ip has recently been identified as the prototype of a new SN IIn class, which is characterized by a history of variability prior to two significant outbursts occurring a few weeks apart (Pastorello et al. 2018). The light curves of SN 2009ip (Margutti et al. 2014) and iPTF13z (Nyholm et al. 2017) also exhibited variability during their main outburst. This might be related to the final explosion of the star, although it is a matter open to debate in the case of SN 2009ip (e.g., Fraser et al. 2013).

The variety of SN IIn properties suggests the presence of multiple progenitor channels. LBVs have been invoked as possible progenitors, as they expel significant amounts of mass during their eruptive phases at velocities similar to those observed in SN IIn spectra (Kiewe et al. 2012; Taddia et al. 2013). In the traditional stellar evolution scenario, however, LBVs are not supposed to be the terminal phase in the lives of massive stars, but they only represent a transitional phase as they evolve and eventually become Wolf Rayet (WR) stars.

Smith & Tombleson (2015) argue that LBVs actually arise from relatively isolated environments, suggesting that they might originate from lower mass stars in binary systems. This result has recently been disputed by Aadland et al. (2018), whose work favors the traditional stellar-evolution scenario for LBVs. Red supergiants with superwinds (Smith et al. 2009) and yellow hypergiants (YHGs, Smith 2014) have also been suggested to produce transients, such as SN 1995N (Fransson et al. 2002) and PTF11iqb (Smith et al. 2015). Electron-capture SNe from super-AGB stars have been suggested to produce SNe-like SN 1994W (Mauerhan et al. 2013) and SN impostors that are similar to intermediate luminosity red transient (ILRT) SN 2008S (Botticella et al. 2009). Wolf-Rayet stars are known to produce CSM interacting SNe, such as the so-called SNe Ibn, which exhibit helium features and lack hydrogen features (Pastorello et al. 2008).

In this paper we focus on the long-lasting Type IIn SN 2013L, whose bright emission was observed for more than 1500 days. This object is also characterized by the presence of multicomponent, asymmetric emission lines, late-time NIR and MIR excess, and very blue early spectral energy distributions (SEDs) dominated by UV emission. SN 2013L was also studied by Andrews et al. (2017, hereafter A17), who provide BgVri photometry, with most of the epochs obtained between +30 d and +192 d after discovery. They also obtained six extra R, V, and i band points at later epochs; the last point was obtained at +1509 d and they analyze early optical and NIR spectra. A17 also provide late-time (up to +1509 d) spectra covering the Hα region.

Here, we present additional data and an analysis of SN 2013L based on observations obtained by the Carnegie Supernova Project (CSP-II; Phillips et al. 2019). This includes optical uBgVri-band photometry extending from +2 d after discovery out to +696 d. Furthermore, we present new, multiple epoch NIR YJH-band photometry, extending from +4 d to +887 d and one epoch of K_s band at +739 d. We also present UV photometry obtained from the Neil Gehrels Swift Observatory, extending from +7 d to +43 d (see also de la Rosa et al. 2016). New optical spectra are also presented as are a number of new NIR spectra that were obtained as part of the CSP-II NIR spectroscopy program (Hsiao et al. 2019). Our aim is to understand the properties of the CSM and the underlying progenitor system of SN 2013L. This comprehensive and high-quality dataset enables us to study SN 2013L in detail, which only a handful of SNe IIn in the literature match (Stritzinger et al. 2012; Fransson et al. 2014).

This paper is structured as follows: in Sect. 2 we provide basic information on SN 2013L and its host galaxy. In Sect. 3 we describe data acquisition and reduction. In Sect. 4, we provide an introduction of the SN light curves and colors. This is followed by the presentation and analysis of the spectra in Sect. 5. SEDs and the bolometric light curve are constructed in Sect. 6. The bolometric properties along with the Hα profiles are modeled in Sect. 7, and the results on SN mass loss and dust emission are provided in Sect. 8. Finally, our conclusions are presented in Sect. 9.

2. SN 2013L and its host galaxy

SN 2013L was discovered by Monard et al. (2013) in ESO 216−39 with an unfiltered apparent magnitude of 15.6. First designated as the possible supernova PSN J11452955−5035531, it was classified by the Public ESO Spectroscopic Survey of Transient Objects (PESSTO; Smartt et al. 2015) as a young SN IIn (Morales Garoffolo et al. 2013). The previous nondetection of SN 2013L (with a limiting unfiltered apparent magnitude of 19) dates 18.957 days prior to discovery, which occurred on 2013 Jan 22.025 UT (JD 2456314.525; Monard et al. 2013). We use days (d) relative to discovery to define the phase of the SN throughout this paper.

The coordinates of SN 2013L are RA(J2000.0) = 11^h45^m29$\stackrel{\text{s}}{.}$55, Dec(J2000.0) = −50°35′53$\stackrel{″}{.}$1. The SN is located 20″ east and 24″ south of the galaxy core. A finding chart of SN 2013L in ESO 216−39 is provided in Fig. 1. In Taddia et al. (2015), we computed the deprojected distance of SN 2013L from its nucleus, which turned out to be 1.09 times the galaxy radius (r₂₅), that is, its location lies in the outskirts of the host. The host-galaxy center is characterized by a N2 (Pettini & Pagel 2004) metallicity of 12+log(O/H) = 8.68 ± 0.18 dex, which when assuming a typical −0.47 dex $r_{25}^{-1}$ metallicity gradient (Pilyugin et al. 2004) implies that SN 2013L was located in a subsolar metallicity environment, that is, 12+log(O/H) = 8.17 ± 0.23 dex. This is typical of long-lasting Type IIn SNe, such as SN 1988Z and SN impostors (Taddia et al. 2015).

Fig. 1. Finding chart for SN 2013L, from a combination of all of the BVgri images obtained with the Swope telescope between JD 2456316.7 and JD 2456448.6.

ESO 216−39 is located at z = 0.01738 (A17; NED¹ lists a slightly lower redshift (0.016992 ± 0.000150), which corresponds to a distance of 72${}_{-5}^{+8.5}$ Mpc (A17). The distance modulus is therefore 34.287${}_{-0.156}^{+0.242}$ mag. According to Schlafly & Finkbeiner (2011), the Milky Way extinction is $A_V^{\text{MW}}=0.367\hspace{0.17em}\text{mag}$ when assuming a Fitzpatrick (1999) reddening law that is characterized by R_V = 3.1. The host-galaxy extinction is neglected in A17, and we assume the same in this work. Moreover, a high-resolution spectrum obtained with the Clay (+MIKE) telescope on +10 d (see Sects. 3 and 5) shows no indication of prevalent Na I D absorption or DIB features. Specifically, the DIB 5780 Å feature is not detectable and the Na I D absorption line has a total equivalent width of 0.14 ± 0.07 Å. This suggests minimal host-galaxy dust reddening. Indeed, according to Fig. 9 of Phillips et al. (2013), this EW value implies $A_V^{\text{host}}\lesssim 0.2\hspace{0.17em}\text{mag}$. In the following, by extinction correction, we mean only Milky-Way extinction correction.

3. Data acquisition and reduction

The CSP-II obtained 49 nights of optical uBgVri photometry and 21 nights of NIR (YJH) photometry. Optical imaging was obtained with the Henrietta Swope 1.0 m telescope (+SITe3 direct CCD camera), whereas NIR imaging was mainly obtained with the du Pont 2.5 m telescope (+RetroCam), while also two nights worth of observations were obtained with the FourStar camera (Persson et al. 2013) on the Magellan Baade 6.5 m telescope. In addition, a single, late K_s-band observation was obtained with FourStar. We also obtained photometric observations with the UltraViolet Optical Telescope (UVOT; Roming et al. 2005) on the Neil Gehrels Swift Observatory (Gehrels et al. 2004) in the UV (uw1, um2, uw2) and in the optical (UBV) over the course of 19 different nights.

Reduction of the optical and NIR data was performed as described in Krisciunas et al. (2017). Given that the SN was still present in our most recent attempt to obtain deep host-galaxy images, we did not perform template subtraction on the science images. However, the background at the location of SN 2013L is quite smooth and low, enabling us to robustly measure its magnitude evolution. Point-spread function (PSF) photometry of the SN was computed and calibrated relative to a local sequence of stars in the field of SN 2013L. The local sequence itself was calibrated relative to the standard-star fields observed over multiple photometric nights from Landolt (1992) (BV) and Smith et al. (2002) (ugri). The NIR J- and H-band local sequences were calibrated relative to the standard stars from Persson et al. (1998), while the Y-band local sequence was calibrated relative to Y-band magnitudes of Persson et al. (1998) standards presented in Krisciunas et al. (2017). Finally, K_s-band photometry of SN 2013L was calibrated relative to 2MASS stars. Photometry of the local sequences in the standard system are listed in Table A.1 (optical) and Table A.2 (NIR), while photometry of SN 2013L in the CSP-II natural system is listed in Table A.3 (optical) and Table A.4 (NIR).

Swift images were reduced using aperture photometry and subtraction of the underlying host-galaxy count rates following the method described by Brown et al. (2014) and by using the zero-points of Breeveld et al. (2011). The Swift UVOT images and photometry are available in the Swift Optical Ultraviolet Supernova Archive (SOUSA; Brown et al. 2014), and they are listed in Table A.5.

The resulting light curves of SN 2013L are shown in Fig. 2. The BgVri photometry and two epochs of MIR photometry that are presented in A17 are also included in Fig. 2. SN 2013L is one of the best observed SNe IIn to date, with data extending from UV to MIR wavelengths.

Fig. 2. CSP-II uBgVriYJHKs-band and Swiftuvw2, uvm2, uvw1-band light curves of SN 2013L (filled symbols). Light curves were shifted for clarity by the amounts shown in the legend. The photometry is also included from A17 (open symbols). Solid lines correspond to tension spline fits to the CSP-II photometry. Epochs of spectral observations are indicated with vertical segments (NIR spectra at the top and optical spectra at the bottom; in black our new spectra, in red those already published by A17). The nondetection limit in an unfiltered/r band is indicated by a red triangle, while unfiltered/r band discovery and confirmation magnitudes are shown by red stars. The top-right inset shows a zoomed in view of the light curves over the first 150 days of evolution.

We used all of the available optical spectra of SN 2013L. This includes four spectra from the ESO Very Large Telescope (VLT) equipped with the X-shooter spectrograph and a single NTT (+EFOSC) spectrum. To these, we added a single spectrum obtained with the du Pont (+B&C spectrograph) telescope and a single medium-resolution spectrum taken with the MagE spectrograph (Marshall et al. 2008) on the Magellan Clay telescope. A single high-resolution spectrum was also taken with the Magellan Clay telescope that is equipped with the Magellan Inamori Kyocera Echelle (MIKE) spectrograph (Bernstein et al. 2003). This spectrum was used to check the host extinction (see Sect. 2) and the early CSM velocity (see Sect. 5). We also used the four late-time spectra presented by A17. In addition to the four X-shooter spectra covering the NIR, we added six previously unpublished NIR spectra obtained with the Magellan Baade telescope equipped with the Folded Port Infrared Echellette (FIRE; Simcoe et al. 2013). Our full NIR spectral time-series extends from +11 d to +861 d.

Optical spectra were reduced in the standard manner, including bias and flat corrections, wavelength calibration to an arc lamp, and flux calibration to a spectroscopic standard star (see Hamuy et al. 2006). We then absolute-flux calibrated the spectra to the CSP-II photometry by performing synthetic (r in the optical and J in the NIR, except the last spectrum where we used the H band) photometry to compute the appropriate scaling of each spectrum in order to match the observed corresponding broad-band magnitude. The X-shooter data were first processed using the esorex pipeline and then extracted and telluric calibrated using our own python scripts. The FIRE spectra were reduced using the firehose software package developed by Simcoe et al. (2013). The reduction steps are described by Hsiao et al. (2019). A log of spectroscopic observations is given in Table A.6 and the full optical and NIR time-series are plotted in Figs. 3 and 4, respectively.

Fig. 3. Visual-wavelength spectroscopic time-series of SN 2013L. Indicated next to each spectrum is the phase, i.e., days from discovery. The spectra that were already published in Andrews et al. (2017) and publicly available are marked by a star next to their phases, the four late-time spectra from A17 are marked by two stars. Each spectrum was normalized by its median and shifted by a constant for better visualization.

Fig. 4. NIR spectroscopic time-series of SN 2013L obtained with FIRE and X-shooter. The phase of each spectrum relative to discovery is indicated next to each spectrum. The spectra that were already published in Andrews et al. (2017) are marked by a star next to their phases. Each spectrum was normalized by its median and shifted by a constant for a better visualization. We marked the telluric bands at the top.

4. Photometry

The definitive photometry of SN 2013L is plotted in Fig. 2. The UV to NIR light curves spanning a total of 13 filters (uw2, um2, uw1, u, B, g, V, r, i, Y, J, H, K_s) are displayed from bottom to top. To appreciate the density of early coverage, the first 200 days of evolution are shown within an inset and also include photometry published by A17. To facilitate a good match between our photometry and that of A17, we introduced small shifts to the latter; specifically, we added 0.1, 0.15, 0.2, 0.12, and 0.2 mag to the BgVri-band light curves, respectively. The shifted light curves agree well with ours at each epoch and in each filter.

The CSP-II light curves commenced within two days after discovery and cover the flux evolution with a high cadence as the SN reached maximum light. The time of maximum occurs first in the bluer bands followed in time with the redder bands. The peak was not observed in the UV bands from uw2 to u due to the rapid drop in UV-emission. The light curves are fit with tension splines to better characterize their shape as well as to infer the epoch and magnitude at maximum light. The epochs of peak and their apparent maximum magnitudes are listed in Table 1.

After the peak, the optical light curves show a rapid decline until about 80 days, when they begin to decline more slowly. The same appears to be true for the NIR, where we have less data covering a 75 day period after the peak. We report the magnitude difference between 80 days and the peak (Δm₈₀) in the optical and NIR in Table 1, as well as the magnitude difference between +80 d and +300 d (Δm₈₀₋₃₀₀). The early decline is faster in the bluer bands. We note that the UV light curves from Swift, which extend to +43 days, also decline very rapidly after the peak and even more rapidly than the u-band light curve. The later decline is different in the different bands; the r band declines more slowly due to the prevalent Hα emission. The g-, r-, and i-band decline after 300 days is faster than the rate measured between +80 d and +300 d. However, after about 700 days when the SN was followed in the R band by A17, the emission appears to flatten out. SN 2013L has a long-lasting light curve, similar to a few other SN IIn events in the literature, and in the following we compare SN 2013L to the data sets of these similar objects.

Using the light curves and their peak apparent magnitudes from Table 1, the extinction, and the distance from Sect. 2, we computed the absolute magnitude light curves and their peaks, which we plotted in Fig. 5 and report in Table 1, respectively. With a peak r-band magnitude of −19.06${}_{-0.16}^{+0.24}$ mag, SN 2013L is among the most luminous long-lasting SNe IIn in the optical, together with SN 1988Z (−19.26 ± 0.3 mag; Turatto et al. 1993) and SN 2010jl (−20.12${}_{-0.19}^{+0.17}$ mag; Fransson et al. 2014; Gall et al. 2014). The early (first ∼130 days) decline rate is rather similar to that of SN 1988Z, SN 2005ip, SN 2006jd (Stritzinger et al. 2012), and iPTF13z (Nyholm et al. 2017), whereas SN 2010jl and KISS15s (Kokubo et al. 2019) decline slower and SN 1996al declines much faster. Later on, in particular after 300 days, the r-band light curve of SN 2013L declines linearly at a similar rate compared to SN 2010jl (which is 2 mag brighter at that phase) and SN 1996al (which is 4 mag fainter at that phase). SN 2006jd and iPTF13z develop a bumpy light curve (only in the r band for SN 2006jd), while the light curves of SN 2005ip indicate nearly constant emission over a period of hundreds of days. During the refereeing process, a work on SN 2015da was published (Tartaglia et al. 2020) which shows another long-lasting SN IIn, 2015da, with peak magnitude and light curve shape similar to those of SN 2010jl.

Fig. 5. Absolute magnitude light curves of SN 2013L covering UV to NIR wavelengths as well as those of the long-lasting Type IIn SN 1988Z (Turatto et al. 1993), SN 1996al (Benetti et al. 2016), SNe 2005ip and 2006jd (Stritzinger et al. 2012), SN 2010jl (Fransson et al. 2014), iPTF13z (Nyholm et al. 2017), and KISS15s (Kokubo et al. 2019).

The NIR (YJH) emission peaks around −19.3 to −19.5 mag, which is rather luminous for a SN IIn; it is 1 mag brighter than SNe 2005ip and 2006jd in the Y and J bands, but fainter than SN 2010jl by about 1 mag. In the H band, SN 2013L declines rather quickly when compared to SNe 2005ip and 2006jd over the first +600 d. In fact, these other long-lasting SNe IIn are characterized by a strong NIR excess due to dust emission from just a few weeks after discovery, which is not immediately evident in SN 2013L over the first +130 d. We return to this point in Sect. 6. In SN 2010jl, the fast optical decline after +400 d corresponds to an increase in the H band, but this is not observed in SN 2013L.

The g − r color of SN 2013L plotted in Fig. 6 evolves from blue to red colors beginning from the onset of our first day of observations and continues to do so to almost +500 d. This evolution is rather fast until +150 d, after which the color evolution slows down. The same is seen in the B − V color curve up to +150 d; however, after +150 d, the B − V color becomes slightly bluer. The g − r color evolution of SN 2013L is similar to that of SN 2005ip, whereas the B − V evolution is similar to that of SN 1996al. On the other hand, SN 2006jd, SN 1988Z, and SN KISS15s show a different evolution, evolving toward the blue over the first +200 d (+400 d for KISS15s). The i − H color evolution of SN 2013L, which is shown in the bottom panel of Fig. 6, confirms that SN 2013L exhibits a smaller NIR excess compared to SN 2005ip and SN 2006jd as well as that SN 2013L is more comparable to the evolution of SN 2010jl at epochs earlier than ∼+450 d. The similarities of its optical-NIR color evolution to objects, such as SN 2010jl, is consistent with dust emission.

Fig. 6. B − V, g − r, and i − H color evolution of SN 2013L as compared to that of other long lasting SNe IIn (SN 1988Z from Turatto et al. 1993, SN 1996al from Benetti et al. 2016, SNe 2005ip and 2006jd from Stritzinger et al. 2012, SN 2010jl from Fransson et al. 2014, iPTF13z from Nyholm et al. 2017, KISS15s from Kokubo et al. 2019).

5. Spectroscopy

The visual-wavelength spectral time series of SN 2013L extending from +6 d to +1509 d is plotted in Fig. 3. The four later spectra, between +482 d and +1509 d, are from A17 and they only cover the wavelength region centered around Hα. The spectra obtained out to +48 d exhibit a blue continuum superposed with narrow Balmer emission lines, which are the only features departing from a smooth black-body (BB) shape. While the spectra become redder, the Hα emission feature grows significantly in strength, emerging by +27 d and more clearly at +33 d. Similar evolution is exhibited by both the Hβ and Hγ emission features. The Ca II NIR triplet emerges in emission from the continuum by +48 d and grows in strength over time as compared to the continuum. In summary, SN 2013L shows typical SN IIn spectral properties at optical wavelengths.

The NIR spectra (see Fig. 4) also cover a large time span, extending from +11 d to +861 d. Until +37 d, they exhibit a BB continuum with narrow Paschen emission features. By +48 d, broad He I λ10830 and Pβλ12822 features emerge. By +144 d, the spectral continuum shows a clear flattening toward the red, which becomes quite prevalent in the last spectrum. These late-time NIR spectra present a number of emission lines in addition to Pβ and He I, which we discuss below.

In Fig. 7 we show a closer inspection of the various spectral features. To do so, we selected spectra after +109 d, where the lines begin to dominate over the continuum and become easier to identify. At optical wavelengths, all the Balmer lines down to Hδ are clearly identified and characterized by broad and similarly shaped emission profiles; see also Figs. 8 and 9. At the blue end of the spectrum, relatively broad features that are compatible with Fe II λλ5018, 5169 are identified. He I λλ5876, 7065 features are also identified and confirmed by the presence of a conspicuous He I λ10830 feature in the NIR. Weak [O I] λ6300 is observed (see also Fig. 10) as well as [O I] λ5577. O I λ7774 and O I λ8446 are only marginally detected, and they are blended with a strong Fe II line and with the Ca II NIR triplet (see also Fig. 11), respectively. The forbidden [Ca II] λ7293,7324 emission lines are also present. The NIR spectra are dominated by Paschen lines, from Pα down to Pη. We note that Pγ is blended with He I λ10830. Bracket lines from Brγ to Brη are also observed to be much weaker than the Paschen series. A bright, broad O I λ11287 line is observed during late epochs, as is evident in Fig. 10.

Fig. 7. Visual (top) and NIR (bottom) spectral line identifications for SN 2013L. It is important to note that the spectra in the top panel were smoothed (dark lines) for presentation purposes. In the optical, we plotted the expected positions for the Fe II lines (red) for a model with Fe II Lyα and Lyβ pumping, from Sigut & Pradhan (2003). These lines can explain the blue excess in the spectrum at late epochs.

Fig. 8. Left panel: Hα, Hβ, and Pβ line profiles compared to one another and plotted in velocity space. Each line has been corrected for extinction and continuum subtracted. Hβ and Pβ lines were scaled by the constants indicated on the left to match the intensity of Hα. The spectral phases in days since discovery are reported on the right. The shape of the line profiles is remarkably similar at all epochs. Top-right panel: Hα luminosity and continuum luminosity in the Hα region as a function of time for SN 2013L and SN 2010jl (Fransson et al. 2014). Bottom-right panel: Hα/Hβ and Hβ/Hγ flux ratio as a function of time for SN 2013L and SN 2010jl.

Fig. 9. Left panel: continuum-subtracted and extinction corrected Hα profiles (black) fitted by a sum of a Lorentzian (blue) and a Gaussian (green) function to reproduce the broad components. The total fit is shown in red. The spectral phases are shown on the right next to each spectrum. The wavelength range around the narrow component was excluded from the fit. Right panel: velocities and luminosities of the Hα components. The centroid and the FWHM of the Gaussian and of the Lorentzian are shown in the top panels, their luminosities and the total luminosity, as well as the velocities of the narrow component are shown in the bottom panels. Here for the narrow component, we considered the velocity of the absorption minimum assuming zero velocity corresponding to the emission peak, such as in A17. We report our results as well as those of A17 to show the close agreement.

Fig. 10. Left panel: normalized spectra in velocity space around the [O I] λ6300 wavelength regions. The original spectrum in black and the smoothed spectrum in red are shown to better guide the eye to follow the line profile evolution. The spectral phases are reported in red next to each spectrum. The emergence of the oxygen line at 6300 Å at late epochs is evident. We note that the line at −3000 km s−1 is Fe II. Right panel: continuum subtracted and extinction corrected spectra in velocity space around the wavelength region of O I λ11827. This line appears at late epochs and it is well fit by a single Gaussian profile of ∼4000 km s−1 FWHM (shown in red). The phases of each spectrum are reported next to them.

Fig. 11. Continuum-subtracted and extinction corrected Ca II triplet profiles (black) fitted by a sum (blue) of three Hα best-fit profiles from Fig. 9 to reproduce the complex line shape. We followed the approach by A17 and confirmed that this is also a good fit for the Ca II triplet at +302 d and +330 d. The spectral phases are shown on the right next to each spectrum. An excess blueward of the best line fit is compatible with the emerging of O I λ8446, but there is also Fe II that can contribute to the flux in that range.

5.1. Balmer lines

Close inspection of the different hydrogen lines reveals remarkably similar features and they evolve in step with time. This is demonstrated in Fig. 8 (left panel) where the Hα, Hβ, and Pβ lines in the +6 d to +330 d spectra are plotted in velocity space and scaled to match their peaks. In doing so, the spectra were first corrected for extinction and then the continuum was removed through the subtraction of a first order polynomial function matched to the continuum. Clearly, the line profiles and their time evolution are nearly identical. A transition from only narrow to a narrow plus a broad component is observed in the +33 d spectrum and significant asymmetric profiles characterize each of the features. In the top-right panel of Fig. 8, we plotted the total luminosity of Hα versus that of the continuum in the Hα region. The line luminosity increases over time, starting from ∼10⁴⁰ erg s⁻¹, and peaks (∼2 × 10⁴¹ erg s⁻¹) at +109 d. At +330 d, it has almost the same value as the peak. Following that, it decreases until it reaches ∼10⁴⁰ erg s⁻¹ at +1509 d. The continuum decreases with time, and the flux from the line already dominates after +68 d.

In the bottom-right subpanel, we show how the values of Hα/Hβ and Hβ/Hγ increase with time after the first +48 d. The behavior of SN 2013L in terms of Hα and continuum luminosity is rather similar to that of SN 2010jl (Fransson et al. 2014), and this also applies when we look at the Balmer line ratios, which we can see in the lower-right panel of Fig. 8.

We further characterize the line profiles to derive information on geometry and velocities in Fig. 9. Following A17, and using Hα after extinction correction and continuum subtraction, we fit the profiles with the sum of two functions to reproduce the broad component. A Lorentzian is fit to the part of the profile symmetric with respect to the zero velocity, and a Gaussian is fit to the bluer part of the profile. The sum of these two functions gives a perfect fit to the profiles, as already shown by A17. We confirm that this also works on our +330 d spectrum, which is not covered by A17. We show the best fit in the left panel of Fig. 9, and the evolution of the velocities that were derived by these fits are in the right panel. The Lorentzian, which is centered around zero velocity at all times, shows a FWHM that increases from ∼2000 km s⁻¹ to ∼5000 km s⁻¹ in the time range between +27 d and +109 d. Then the FWHM velocity of the Lorentzian component drops to ∼3000 km s⁻¹ at +330 d and down to ∼1000 km s⁻¹ at +1509 d. The Gaussian, reproducing the blue part of the Hα profile, is blue-shifted by ∼4000 km s⁻¹ at early epochs. The blue-shift slowly decreases down to ∼1000 km s⁻¹ at +1509 d. The FWHM of the Gaussian drops with time from ∼6000 km s⁻¹ at early epochs to ∼2500 km s⁻¹ by +330 d and later epochs.

The Hα total luminosity, which we describe above, is dominated by the broad Lorentzian component, and both the Lorentzian and the Gaussian luminosity behave similarly. The narrow component of Hα (see Fig. 9, bottom-right panel) has a decreasing minimum absorption velocity ($V_{\min }^{\mathrm{P}-\mathrm{Cyg}}$, measured by A17 and by us, assuming a velocity of zero at the emission peak of Hα) between −68 km s⁻¹ and −130 km s⁻¹.

For an estimate of the CSM velocity from the narrow Hα component (V_w), one should use the blue-velocity-at-zero-intensity (BVZI) of this narrow component, assuming zero velocity at a Hα-rest wavelength, as in Fig. 12. This velocity is constant around 120 km s⁻¹ during the first +109 d. It becomes higher (up to 240 km s⁻¹) at later epochs. The resolution of the spectra (see A17, their Table 3) is sufficient to resolve these narrow lines to obtain their velocities (FWHM resolutions of X-shooter and IMACS are 40 and 57 km s⁻¹, respectively). Values of wind velocity in the range of 120−240 km s⁻¹ are also observed in other SNe IIn (see, e.g., Stritzinger et al. 2012; Taddia et al. 2013; Fransson et al. 2014; Gall et al. 2014).

Fig. 12. Narrow component of Hα plotted in velocity space. The spectral phases in days since discovery are reported on the right. The zero velocity was chosen to correspond to 6562.8 Å in the rest frame. We plotted (in red) the best fit to the line profile, which was achieved by fitting two Gaussian components, one for the absorption and one for the emission. The emission dominates at early epochs, while absorption becomes stronger at later epochs. We marked the BVZI with a dashed segmented line that ranges from −120 to −240 km s−1.

The first spectrum where we were able to measure the narrow P Cygni minimum is the high-resolution MIKE spectrum taken at phase +10 d. The NTT (+6 d) and du Pont (+27 d) spectra do not have the resolution to measure the P Cygni absorption profile. Figure 12 also shows a transition in the shape of the narrow Hα, with the emission part dominating at early phases and decreasing with time until the absorption component dominates after +330 d.

The decomposition of the broad profile of Hα into one Lorentzian and one Gaussian component reproduces the line shape well, which was already shown by A17 and is confirmed here. However, we also applied a different, more physically-motivated decomposition to the line. This is based on the observation that the Hα profile can also be the sum of a broad, flat, skewed component, which decreases its FWHM with time, and of a narrower component centered at zero velocity whose relative strength increases over time. This is well illustrated in Fig. 13, where we normalized the continuum-subtracted Hα profiles to the flux of the “shoulder” of the broad, skewed component. This shoulder moves toward lower velocities with time. This decomposition is suggested above all by the spectral modeling we present in Sect. 7.2, which reproduces the broad, flat component with the emission from a narrow emitting shell. While this would give a boxy profile if the medium is optically thin, electron scattering smooths this profile and occultation results in a skewed blue-shifted profile (Sect. 7.2).

Fig. 13. Evolution of the continuum subtracted Hα line profile in velocity space. The flux was normalized to the flux of the flat part of the profile between the shoulder and the central component. The spectra from +482 d to +1133 d are from A17.

In Fig. 14 (left panel) we fit the Hα profiles as a combination of a broad skewed Gaussian with two different HWHM (half-width-at-half-maximum) (blue) and a central, narrower Gaussian (green) function. The central Gaussian in the first spectrum is too weak to be seen. Here, we ignored the narrow P Cygni region, as in Fig. 9. The centroid of the two components and their FWHM and HWHM are shown in the top panels. The narrower Gaussian is basically centered at zero velocity, while the skewed Gaussian moves to lower velocities from +100 d to +1500 d. The FWHM of the broad component decreases with time, except for the first two epochs where it is more difficult to measure. The FWHM of the central Gaussian is rather constant, but it decreases in the last two spectra. The luminosities of the two components behave similarly, rising to a peak and then declining. In the bottom-right panel, we report the BVZI of the narrow P Cygni, which is discussed above and shown in Fig. 12.

Fig. 14. Continuum-subtracted and extinction corrected Hα profiles (black), this time fitted by a sum of a skewed, broad Gaussian (blue) and a narrower Gaussian (green) function to reproduce the broad components. The total fit is shown in red. The spectral phases are shown on the right next to each spectrum. The wavelength range around the narrow component was excluded from the fit. Right panel: velocities and luminosities of the Hα components. The centroid and the FWHM and HWHM of the skewed broad Gaussian (blue) and of the central Gaussian (green) are shown in the top panels, their luminosities and the total luminosity, as well as the velocity of the narrow component are shown in the bottom panels. Here, for the narrow component, we consider the blue-velocity-at-zero-intensity (BVZI) assuming zero velocity at the Hα rest wavelength to obtain a better estimate of the CSM velocity.

Finally, we measured the fluxes of three bright hydrogen lines from the NIR spectra and we show their similar evolution in Fig. 15. These lines rise rapidly in flux to reach their peak between +100 d and +200 d.

Fig. 15. Flux as a function of time for three bright hydrogen lines in the NIR spectra. Their evolution is rather similar.

5.2. Lines from heavier elements and continuum

Other interesting features present in the spectra of SN 2013L are the oxygen lines. In Fig. 10, we see broad [O I] λ6300 beginning from +109 d and lasting onward. The same is true for O I λ11287, which becomes very strong by +302 d. The FWHM velocity of O I λ11287 at this phase is around 4300 km s⁻¹, which is obtained by fitting a Gaussian to the continuum-subtracted line after extinction correction and absolute calibration. The O I λ11287 persists until the last spectrum at +861 d. The FWHM of the O I λ11287 line almost matches the FWHM of the Hα profile (see Fig. 16).

Fig. 16. Continuum-subtracted Hα and O I λ11287 at +300 to +330 d. The NIR O I line, which was scaled to match the peak of Hα, shows a broad profile that is almost as broad as Hα.

Following A17, the Ca II NIR triplet from +68 d to +330 d (shown in Fig. 11) is well-represented by a combination of three Hα profiles (one for each of the three Ca II lines), which were scaled by three different constants to fit the line profile. In the spectra at +302 d and +330 d, there is a small excess in the bluer part of the line that might be compatible with O I λ8446, which emerges at the same time as O I λ11287 and also O I λ7774 (see Fig. 7). It might also be compatible with a Fe II line, as shown in Fig. 7. Also the blue excess between 4000 Å and 5400 Å can be explained by a forest of Fe II lines (see Fig. 7). Neither of the O I or Ca II lines show any velocity shift and have a FWHM that is compatible with the central component of Hα. They are therefore likely to arise in the CSM.

In addition to the line analysis, it is of interest to study the evolution of the spectral continuum, especially in the NIR. Plotted in the inset of Fig. 17 are all the optical spectra after applying an extinction correction and an absolute flux calibration relative to r-band photometry; they were then fit with a BB function. Here we exclude the region of the spectra dominated by strong emission lines. A single BB function provides a good representation of the continuum. In the NIR, where we observed spectra to later epochs, a single BB function provides a reasonable fit until +109 d. However, the +144 d spectrum shows a better fit for a two-component BB function, with a hot component fitting the bluer portion of the spectrum and a warm component dominating the flux above 20 000 Å. The same is true for the later spectra, where the warm BB component becomes relatively stronger. In the last NIR spectrum, the second BB component fitting the red portion of the spectrum is so strong that the continuum is more prevalent above 20 000 Å than in the blue. From the BB fit, estimates of the temperature, radius, and luminosity were derived (see Sect. 6) together with the same quantities derived from fitting SEDs constructed from broad-band photometry.

Fig. 17. Inset: optical spectra after absolute flux calibration and extinction correction, fitted by a single BB function. The parts of the spectrum in light gray were excluded from the fit because of the presence of strong emission lines. Main panel: NIR spectra after absolute flux calibration (using J band) and extinction correction, fitted either by a single BB (at early epochs), or by a sum (magenta) of two BB components, one warm (red) and one hot (blue). The warm (red) component dominates at late epochs. The flux of the last spectrum was multiplied by five for a better visualization.

6. SEDs and constructing the bolometric light curve

The spectral continuum fit discussed in Sect. 5.2 and shown in Fig. 17 reveals an interesting double BB component emerging at late epochs. We now turn to photometry and build SEDs extending from UV to NIR wavelengths. From these, we further examine the different underlying emission regions and also construct a bolometric light curve.

Using the tension splines fitted to the photometry presented in Fig. 2, we obtained interpolated UV to NIR photometry from +7 d to +335 d, when photometry was obtained in all of the CSP-II bands extending from u to H. Between +7 d and +43 d, we also obtained interpolated UVOT photometry. Next, appropriate AB offsets that transform the CSP-II natural system to the standard system were added to the photometry (see Krisciunas et al. 2017, their Table 16). The photometry was then corrected for reddening and then converted into monochromatic fluxes at the effective wavelength of each filter. Six examples of our SEDs are plotted from +10 d to +335 d in Fig. 18. Here, we selected epochs where the NIR fluxes are directly available from the observations, and not just interpolated. Up through +68 d, the SEDs clearly exhibit a shape that is well-matched by a BB function, with the exception of the r band (red diamonds), which begins to deviate from the best BB fit due to prevalent Hα emission. The r-band flux was therefore not included in the BB fits, which are shown with a black line. The UVOT fluxes in the UV (white, black, and gray symbols) nicely follow the BB shape, indicating that at least at early epochs, the UV is not strongly suppressed by line blanketing. The NIR fluxes (light gray, green-cyan, and cyan points) fall on the BB Rayleigh-Jeans tail until +68 d; however, at later epochs such as +335 d, they show an excess that is reproduced well by a second BB function (red line). The BB temperature drops from an initial 10 700 K at +10 d to 5800−6100 K between +68 d and +335 d. The second component in the NIR has a lower temperature of about 1400 K at +335 d. At +335 d, the r-band flux point is well above the SED continuum as a result of the prevalent Hα feature. We performed BB fits to the SEDs from +7 d to +335 d, and we notice that the second BB component gives a better fit to the SED only starting from +132 d, which is in agreement with the spectral fits shown in Fig. 17, which do not show a second component at +109 d but do so at +144 d.

Fig. 18. Example of SEDs of SN 2013L with the specific epochs corresponding to times when NIR photometry is available. Each flux point was obtained by converting the magnitudes into fluxes after having applied the appropriate extinction correction. The color and symbol code for each band is the same as in Fig. 2. A single BB function was fit to each SED until +68 d (black line), while a sum of two BB functions was used to fit the +335 d SED in order to reproduce the NIR excess (red line) and the optical BB shape (blue line). The r band is not included in the fit due to the prevalent Hα emission. The BB temperature is reported in each subpanel as well as the phase of the SED.

At very late epochs (+630 d to +850 d), when only K_s band and Spitzer photometry at 3.6 and 4.5 μm (see A17) are available, the NIR (and MIR) excess is even stronger, as already shown in the last NIR spectrum in Fig. 17. This excess is highlighted in Fig. 19, where the full SED is shown, including all of the flux points and our last NIR spectrum. A sum of two BBs is confirmed to match the emission including the MIR flux points.

Fig. 19. Late phase Ks-band and Spitzer photometry converted into monochromatic fluxes plotted together with our last epoch NIR spectrum. The SED exhibits a strong NIR and MIR excess that dominates the SN emission two years after discovery. The NIR spectrum was scaled to match the Ks-band flux. We fit a sum of two BB functions (magenta dashed line) to reproduce the NIR spectrum and the MIR fluxes from the first Spitzer epoch. The warm BB component (red dashed line) dominating the NIR and MIR has a temperature of 1022 K and the hot component (blue dashed line) has a temperature of 4513 K.

The BB fits to the SEDs are useful to build a bolometric light curve. In order to obtain the luminosity, the BB fluxes were multiplied by 4πD² where D is the SN distance. The luminosity from the two BB components is shown in the top panel of Fig. 20, as well as their sum. The warm (NIR) component emerges only after +132 d. We also show the luminosity when adding the remaining Hα flux as the r-band monochromatic flux was excluded from the BB fit. For comparison purposes, we show the luminosity derived from the BB fit to the spectra. The optical spectral fits give very good agreement in luminosity, the NIR spectral fits are not accurate in reproducing the hot BB luminosity in the optical given the different wavelength coverage, but they provide an estimate for the warm BB luminosity in the NIR.

Fig. 20. Top panel: bolometric luminosity of SN 2013L. From the BB fit to the photometry (solid lines), we derived the luminosity of the hot BB component (solid blue line, dominating at early epochs) and of the warm BB component (red solid line, dominating at late epochs and emerging at +132 d). Their sum (dashed magenta line) and their sum plus the Hα emission neglected from the BB fit (black solid line) are also shown, the last one represents the total luminosity. We also obtained the corresponding luminosities from the BB fit to the optical (blue diamonds) and NIR spectra (blue dots and red squares). Central panel: temperature evolution of SN 2013L for the two different BB components fitting the SEDs of SN 2013L (symbols as in the top panel). In addition, the dashed lines mark the uncertainties on the temperature from the BB fit to the photometry. Bottom panel: radius evolution for the two different BB components fitting the SEDs of SN 2013L, symbols are the same as those in the top and central panels.

Even though it was discovered before optical peak, due to the huge UV flux, the bolometric luminosity at early epochs indicates that SN 2013L was discovered after the bolometric peak. The observed bolometric peak is at ≳3 × 10⁴³ erg s⁻¹. The luminosity drops to 2 × 10⁴² erg s⁻¹, one order of magnitude, in ∼100 days. Subsequently its luminosity declines more slowly. The warm BB luminosity contributing to the total luminosity is rather constant with time, at about 1 × 10⁴² erg s⁻¹, and it begins to dominate at roughly +300 d to +400 d after discovery. Then it largely dominates by +800 d when the hot, optical BB luminosity is 3−4 × 10⁴⁰ erg s⁻¹ and the warm, NIR BB luminosity is at 7 × 10⁴¹ erg s⁻¹.

Plotted in the central and bottom subpanels of Fig. 20 are the temperatures of the hot and warm BB components obtained from our BB function fits to our spectrophotometry. The hot BB temperature drops from 15 000 K right after discovery to about 6000 K at +1000 d, and then it stays constant. The NIR warm BB temperature is at ∼1000−1700 K from +132 d until +861 d.

The BB radius for the hot component expands until +70 days (from 1 × 10¹⁵ cm to 2.5 × 10¹⁵ cm), then it begins to slowly decrease, reaching 3 × 10¹⁴ cm at +800 d. The warm NIR BB radius stays rather constant at 1−3 × 10¹⁶ cm.

Although the BB fits give a good agreement with the observed SED, this may, especially at later epochs, underestimate the luminosity and also the radius of the hot component. This is illustrated by the simulation of SN 2010jl in Dessart et al. (2015, their Fig. 13), where the model has a considerable UV contribution compared to a BB fit of only the optical range. This in turn leads to a large underestimate of the bolometric luminosity. We return to this in Sect. 7.4. In addition, when scattering dominates, the radius of the photosphere is significantly different from that derived from a BB-fit.

7. Modeling

In the following, we detail a model consisting of various emission features driven from the interaction between the rapidly expanding supernova ejecta with circumstellar material. One of the goals of the ascribed model is to estimate the mass-loss rate and hence to clarify the progenitor scenario. To begin, the Hα profile was modeled; in the following, a model of the bolometric light curve is provided, making use of information derived from the spectra.

7.1. Explosion epoch

In order to accurately apply a CSM-interacting model to describe SN 2013L, it is first required to estimate its explosion epoch. We know that the last nondetection occurred 19 days before discovery (see Sect. 2) and given the magnitude limit (19 mag), the SN certainly exploded within this 19 days interval. Considering the radius provided in Fig. 20, we can extrapolate to r_BB = 0 by fitting the first +10 d after discovery with a straight line since its evolution appears linear. This provides an estimate of the explosion epoch of −15.0 d prior to the discovery epoch, which is within the interval determined by the last nondetection epoch. If instead of fitting only up to +10 d, we rather fit to +30 d with a second order polynomial (the evolution resembles a parabola), we obtain an explosion epoch on −15.6 d prior to discovery, that is, very similar to the first estimate. In the following, the explosion epoch is inferred to have occurred on −15 d, that is, JD_explo = 2456299.5${}_{-4}^{+15}$. The radius fit is shown in Fig. 21. A caveat to our explosion time estimate is the possibility of a large (10¹⁴−10¹⁵ cm) progenitor radius, which would imply a later explosion date.

Fig. 21. Fit of the radius evolution used to derive the explosion epoch. A linear fit from the first 10 days (red) gives an explosion epoch of −15 d before explosion. A similar explosion epoch was obtained by fitting the first +30 d with a second-order polynomial (magenta).

7.2. Modeling of the spectral lines

The main characteristics of the spectra are the dominance of the Balmer lines, the wide line profiles, and the blue-shift of the lines. To illustrate the evolution of the line profiles, we show the Hα line from +48 d to +1133 d after explosion in Fig. 13. To clearly see the evolution of the shoulder, we normalized the flux to that of the flat part of the shoulder.

Figure 13 clearly demonstrates the gradual decrease of the velocity of the shoulder in the blue wing of the line, from ∼4000 km s⁻¹ at early epochs to a near disappearance of the shoulder at the last epochs. The figure also shows the separation of the line profile via two distinct components, a flat part with decreasing blue-shift and a central nearly symmetric line, with a fairly constant FWHM.

To model the dynamic line profile, we note that a flat-topped line profile is indicative of a hollow emitting shell, moving with a high expansion velocity. The observed profiles are, however, asymmetric and smooth with very extended wings, which are characteristic of Thomson electron scattering. To model the line profile, we used a Monte Carlo code, calculating the emitted spectrum after scattering for a given input spectrum. The code is a modified version of the code used for modeling the line profiles of SN 2010jl (Fransson et al. 2014); it calculates the emission including the scattering by the thermal electrons as well as a macroscopic bulk velocity. We emphasize that we do not aim to develop a complete self-consistent model here. This requires a full radiation-hydrodynamical code, which is beyond the scope of this paper. Here we mainly aim to use the line profiles and light curve as diagnostics of the main properties of the SN and, in particular, the shock velocity and mass-loss rate.

We assume that the Hα photons are emitted from a spherical shell with an inner radius of R_em in and an outer radius of R_em out, and that the velocity inside R_em out is homologous, that is, V(r) = r/t. For example, a constant velocity shell would give a similar result. The natural interpretation is that the emitting shell is the cool, dense shell (CDS) resulting from the interaction region between the expanding ejecta and the dense CSM, as discussed in Chevalier & Fransson (1994) or Dessart et al. (2015) for SN 2010jl. The photons are then followed as they are scattered by the slow, ionized, but unshocked CSM outside of the shock front. In the standard scenario, we assume this CSM to have a density ∝r⁻², which is characteristic of a constant velocity wind before the explosion, but we also discuss cases with a different dependence, which may apply for a time limited eruption with a varying mass-loss rate. Because the shell is narrow, the latter assumption is not very important. We discuss the density and temperature structure further below.

The ionized CSM is mainly characterized by an optical depth τ_e and a temperature T_e. Because the dispersion in frequency by the incoherent electron scattering is $\mathrm{\Delta }V\approx N_{\mathrm{scatt}}^{1/2}{V}_{\mathrm{thermal}}\propto {\tau }_{\mathrm{e}}{T}_{\mathrm{e}}^{1/2}$ (since $N_{\text{scatt}}\approx {\tau }_{\text{e}}^2$), the optical depth and temperature of the CSM are not independent. A different average CSM temperature would therefore give a similar result for the line profile for an optical depth scaling as τ_e ∝ (T_e/10⁴ K)^−1/2.

In realistic models, the temperature depends on the distance from the shock. The calculations in Dessart et al. (2015) show a temperature decreasing from ∼2 × 10⁴ K to ∼1 × 10⁴ K for the SN 2010jl model at 30 days. Here, we model this by a temperature decreasing as T_e ∝ r^−0.4 from the shock, but the line profile is not very sensitive to this assumption. Also for other parameters, such as the bound-free opacity, we use the model by Dessart et al. (2015) as a guide. It is important to note, however, that this model was used for a relatively low ejecta mass of 9.8 M_⊙, while a more massive ejecta may give different parameters.

The importance of scattering in the post-shock region and ejecta depends on the degree of thermalization there. Because of the high density and low temperature in the cooling gas in the CDS, the bound-free and free-free continuum absorption increases compared to electron scattering, as can be seen in the simulations of Dessart et al. (2015). The exact level depends on the shock compression, which is sensitive to multidimensional effects, as well as magnetic fields and the specific ejecta model. The details of the shock properties and the radiation from this are, however, not yet fully understood (see e.g., discussion in Waxman & Katz 2017).

In assuming a viscous shock and ignoring preacceleration of the CSM, the temperature of the shock is

$$\begin{array}{c}\hfill T_{\mathrm{s}}=3.4\times {10}^8{\left(\frac{V}{5000{\mathrm{km}\mathrm{s}}^{-1}}\right)}^2\mathrm{K},\end{array}$$(1)

for solar abundances. Dessart et al. (2015), however, find a considerably lower shock temperature of ≲10⁵ K, including the effects of radiation, although they note that the shock structure is not fully resolved in their standard calculation, while higher temperatures are obtained in high resolution simulations. The presence of hot gas with T_e ≳ 10⁸ K from a viscous shock is consistent with the X-ray observations of SN 2010jl (Ofek et al. 2014b; Chandra et al. 2012, 2015).

For a CSM with a constant mass-loss rate, Ṁ, and with a wind velocity, u_wind, the density corresponds to ρ(r) = Ṁ/(4πu_windr²). For a more general CSM with a power law density ∝r^−s, this corresponds to

$$\begin{array}{c}\hfill \rho (r)=\frac{\dot{M}}{4\pi u_{\mathrm{wind}}{R}_0^2}{\left(\frac{R_0}{r}\right)}^{\mathrm{s}},\end{array}$$(2)

where r₀ is a reference radius. The radius of the shock is then

$$\begin{array}{c}\hfill R_{\mathrm{s}}=R_0{(t/t_0)}^{(n-3)/(n-s)},\end{array}$$(3)

assuming a similarity solution can be applied (Chevalier 1982). The shock velocity evolves as

$$\begin{array}{c}\hfill V(t)=V_0{\left(\frac{t}{t_0}\right)}^{(s-3)/(n-s)}.\end{array}$$(4)

Here V₀ is the velocity at a reference time t₀. For a CSM shell with an outer radius R_out, we can take R₀ = R_out, or in terms of the time of the breakout of the shell, t_b, R_out = V(t_b) t_b.

The cooling time of the shock is

$$\begin{array}{cc}\hfill t_{\mathrm{cool}}=& 5.7{\left(\frac{\dot{M}}{0.1M_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-1}\left(\frac{u_{\mathrm{wind}}}{100{\mathrm{km}\mathrm{s}}^{-1}}\right)\hfill \\ \hfill & \times {\left(\frac{V}{5000{\mathrm{km}\mathrm{s}}^{-1}}\right)}^3{\left(\frac{t}{100\mathrm{days}}\right)}^2\mathrm{days},\hfill \end{array}$$(5)

assuming free-free cooling dominates. Compton cooling may decrease this further. The gas therefore rapidly cools until the cooling is balanced by the radiative heating. Because the post-shock gas is in near pressure balance, the compression is ∼T_s/T_ps ≳ 10⁴, where T_ps is the post-shock temperature, set by the radiation field. In reality, instabilities and magnetic and radiation pressure probably limit the compression (Steinberg & Metzger 2018). In any case, the region is cool and dense, and thermalization by bound-free absorption is very important.

The total electron scattering depth of the CSM (if completely ionized, see below) is then

$$\begin{array}{c}\hfill {\tau }_{\mathrm{e}\mathrm{CSM}}=\frac{\dot{M}\kappa }{4\pi (s-1)u_{\mathrm{wind}}V(t_{\mathrm{b}})t_{\mathrm{b}}}[{\left(\frac{t_{\mathrm{b}}}{t}\right)}^{(s-1)(n-3)/(n-s)}-1].\end{array}$$(6)

Here, we ignore the contribution from the lower density CSM outside R_out, which provides a small extra contribution.

The swept up shocked CSM mass is

$$\begin{array}{c}\hfill M_{\mathrm{CDS}}=\frac{\dot{M}}{(3-s)u_{\mathrm{w}}}{\left(\frac{R_{\mathrm{b}}}{R_{\mathrm{s}}}\right)}^{s-2}{R}_{\mathrm{s}}[1-{\left(\frac{R_{\mathrm{in}}}{R_{\mathrm{s}}}\right)}^{3-s}],\end{array}$$(7)

where R_in is the inner radius of the CSM. Also, in assuming this to be fully ionized, the electron scattering depth is

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{e}\mathrm{CDS}}& =\frac{\dot{M}\kappa }{4\pi (3-s)u_{\mathrm{w}}{R}_{\mathrm{s}}}{\left(\frac{R_{\mathrm{b}}}{R_{\mathrm{s}}}\right)}^{s-2}[1-{\left(\frac{R_{\mathrm{in}}}{R_{\mathrm{s}}}\right)}^{3-s}]\hfill \\ \hfill & \approx \frac{\dot{M}\kappa }{4\pi (3-s)u_{\mathrm{w}}{t}_{\mathrm{b}}V(t_{\mathrm{b}})}{\left(\frac{t_{\mathrm{b}}}{t}\right)}^{(s-1)(n-3)/(n-s)},\hfill \end{array}$$(8)

if R_in ≪ R_s.

We treat the structure of the shocked gas as a region of hot gas with T_e = T_s and thickness t_coolV_s/4, and inside of this are a dense shell with optical depth τ_e CDS and a temperature ∼10⁴ K. The optical depth of the hot gas immediately behind the shock is ≲0.1 for V_s ≲ 5000 km s⁻¹ and it is not important for the scattering of the optical photons. Because of the need for thermalization in the CDS, the structure of the ejecta interior to the CDS is not important, and we assume the density to be constant and similar to that in Dessart et al. (2015).

The main free CSM parameters are therefore Ṁ, V(t_b),T_e, and ϵ. Finally, the line profile depends on the Hα emissivity from the shocked gas as well as the CSM. If recombination dominates, we expect the emissivity, j, to scale as $j\propto {\rho }^2{T}_{\text{e}}^{-0.942-0.031\hspace{0.17em}\log (T_{\text{e}}/{10}^4\hspace{0.17em}\text{K})}$ (Draine 2011). At high Hα optical depth and electron density this is likely to be a rough approximation at best, especially in the CDS. We use this expression for the CSM, where the radial dependence of the emissivity is of some importance. However, the small radial extent of the CDS and its dominance for the luminosity means that the radial dependence is unimportant and the main parameter is the total luminosity relative to that of the CSM luminosity. We illustrate the effect of this below.

The shape of the line profile depends on Ṁ, V(t_b), T_e, and ϵ in different ways. The most important aspect for the asymmetry of the line profile and also the main result of this modeling is the velocity of the shock, which determines the velocity shift of the shoulder in the blue wing of the profile.

The optical depth of the CSM, τ_e CSM, is the most important parameter for the slope of the blue wing shortward of the shoulder. Together with V_s, this determines the mass-loss rate needed. There is, however, a degeneracy of the mass-loss rate and the outer radius, or equivalently t_b, for a given V_s (Eq. (6)).

In our line profile calculations, we treated the ratio of absorption to electron scattering as a parameter, ϵ. In the CSM Dessart et al. (2015) find ϵ ∼ 10⁻³; however, in the dense post-shock gas, ϵ depends on the density and temperature in the CDS and it may be ≳1, resulting in strong thermalization of the Hα photons. Thermalization in the CDS and ejecta are important, while for ϵ ∼ 10⁻³ they are less important in the CSM. An argument for high values of ϵ in the CDS and ejecta come from the fact that multiple scattering between regions with high relative velocities, from either the CSM to the ejecta or within the ejecta, produces a strong red wing extending to high velocities (e.g., Fransson & Chevalier 1989). This is not seen in SN 2013L. To suppress this wing, an efficient thermalization in the ejecta is required. In our models, we have assumed ϵ ≳ 1 in the CDS.

In the left panel of Fig. 22, we show the result of a simulation of the observed Hα line at +48 d. The parameters for this simulation are V_s = 4800 km s⁻¹, τ_e,  CSM = 12.3, and τ_e,  CDS = 8.6.

In the right panel, we also show a model for the line at +330 d. The corresponding parameters at this epoch are V_s = 2700 km s⁻¹ and τ_e,  CSM = 7.2. In both simulations, we used T_e = 1.5 × 10⁴ K at the shock, decreasing to 9000 K at the outer boundary of the CSM shell. For the observed line profile, we subtracted the continuum by fitting a straight line between ±15 000 km s⁻¹ from the line center. We did not attempt to model the narrow P Cygni profile of the outer CSM, with velocity ∼120 km s⁻¹.

From these figures, it is seen that this simple model with a few parameters can provide a very good fit to the observed line profile outside of ±3000 km s⁻¹, as shown by the cyan colored lines, which show the difference between the model and the observations. As illustrated by the magenta, dashed lines, showing the input emissivity, the line profile at these dates is dominated by the Doppler broadening from the macroscopic velocity of the CDS. However, in the center of the line, below velocities ±3000 km s⁻¹, there is a significant difference. This residual is interesting since it has a similar width as the central component, which becomes increasingly dominant at the later phases, and is very prominent by +330 d and later. We return to this component below.

The agreement between the model and the observed line profiles shows that an expanding Hα emitting shell, with electron scattering in a dense CSM outside the shell, provides a reasonable fit to the spectrum. While several of the parameters, in particular the temperature, emissivity, and thermalization parameter, ϵ, are model dependent and require a more sophisticated NLTE calculation, as in Dessart et al. (2015), the optical depth of the CSM and, most important for this paper, the shock velocity should be robust. The shoulder in the blue wing of the line therefore gives a measure of the velocity and its evolution with time.

Fig. 22. Continuum subtracted Hα line (shown in black) from days 48 (left) and 330 (right) together with the electron scattering profiles from the model with Hα emission only from the ejecta and shocked region (red). The dashed magenta lines show the box-like input emission profiles from a radially thin shell with velocity 4800 km s−1 and 2700 km s−1, respectively. The cyan lines show the residual between the model and the observed profile. Electron scattering and occultation by the photosphere result in the smooth, blue-shifted profile. We note the central excess from the models, indicating a different component.

In Fig. 23, we show the observed velocity of the shoulder of the line profile for the different dates. For a number of selected dates, including the ones shown in Fig. 22, we also plotted the velocity of the Hα emitting shell as blue squares. For both sets of data, we fit a power law to the velocity as a function of time. For the velocity of the modeled shell before electron scattering, we find that

$$\begin{array}{c}\hfill V_{\mathrm{shell}}=4000{\left(\frac{t}{100\mathrm{days}}\right)}^{-0.23}{\mathrm{km}\mathrm{s}}^{-1}.\end{array}$$(9)

Fig. 23. Red squares: velocity evolution of the Hα-line blue shoulder as function of time. Blue squares: velocity of the shell from the electron scattering model for selected dates. The dashed lines show a least square fit to these.

The slope of the velocity of the blue shoulder has the same time dependence, but a lower normalization, V_break = 3100 (t/100 days)^−0.23 km s⁻¹. This difference is also illustrated Fig. 22, where it is seen that the electron scattering smooths the intrinsic line profile, causing the blue shoulder in the line profile to move to a lower velocity.

7.3. The location of spectral line formation

From the above modeling, it is clear that most of the Hα emission must originate from gas behind the shock. For a dense CSM, this is from the radiative forward shock (e.g., Chevalier & Fransson 2017). This explains the high velocity needed for the majority of the photons, which are scattered by gas in this region and by gas that is external to the shock. The X-ray and UV emission from the shock is also the source of the ionization of the dense unshocked CSM that is needed for most of the electron scattering. This explains the major part of the line profile shown in Fig. 22.

The presence of a dense slow moving CSM is apparent from the narrow P Cygni lines seen on top of the Balmer lines. From the spectra, Hα was measured to have a BVZI of 120−240 km s⁻¹. This line component, as in other SNe IIn, is formed from unshocked CSM that is excited by the emission coming from the shock region. As a result, the BVZI value of this narrow component provides a measure of the CSM velocity. This component is, however, likely to be only the outer region of the CSM, or a separate extended component with lower density, since the P Cygni line would otherwise be smeared out by electron scattering (see below).

The origin of the low-velocity component below ±3000 km s⁻¹ (the central component of Hα), which is apparent as the residual in the simulations of the line profiles in Fig. 22, is especially interesting. The fact that the luminosity of this component approximately follows that of the broad component (see Fig. 14) suggests that a similar excitation mechanism is occurring for both, that is, excitation by the radiation from the shock region. Excitation by radioactivity in the core is very unlikely, both because of the high luminosity, the time evolution, and the fact that the ejecta are opaque to electron scattering and bound-free absorption.

We believe that the most likely explanation is emission coming not only from the high velocity shock, but also from Hα emitting gas ahead of the shock, with low velocity that is characteristic of the unshocked CSM. The motivation for this comes from the symmetric residual in the modeling shown in Fig. 22, which is characteristic in shape of an electron scattering line from a stationary medium. Referring to shock models, this is also a natural consequence of preionization of the unshocked gas by the UV and X-ray emission from the shock, which heats the gas to 10 000−20 000 K (e.g., Dessart et al. 2015). This is partly the same region that is giving rise to the scattering of the Hα emission from the shocked gas, discussed above. Because the emission is $\propto n_{\text{e}}^2$, while the electron scattering depth is ∝n_e, it is the inner denser regions close to the shock that mainly contribute to the emission, while the scattering is done over a larger radius. The main difference between the two emitting components is the low velocity of the unshocked emitting gas. To model this quantitatively, we added an emitting component coinciding with the inner region of the ionized CSM, immediately outside the shock, with emissivity as discussed in Sect. 7.2. The main requirement is that the Hα emitting region should coincide with the electron scattering region and have a low velocity. These requirements are met by the preionization region in front of the shock.

In Fig. 24, we show the result of a model where we have adjusted the ratio of the Hα luminosity of the unshocked and shocked gas to give the observed ratio of the narrow and broad components of the line. The velocity of the unshocked gas has been assumed to be similar to that inferred for the narrow P Cygni absorption, ∼100 km s⁻¹. As is apparent from the figure, this model greatly improves the agreement between the observations and model and now includes the full line profile, with the exception of the narrow P Cygni profile below ∼100 km s⁻¹. The latter must come from gas, which is optically thin to electron scattering, and most likely from gas located outside of the high-density CSM.

Fig. 24. Same as Fig. 22, but now also including Hα emission from the unshocked medium preionized by the shock, which can be seen as the narrow magenta line on top of the broad from the shocked gas (normalized to one at maximum). The magenta, dashed profile shows the line profile without electron scattering or continuum thermalization. We note that the residual close to the line center is now much lower at both epochs. The steep profile on the red side is a result of the comparatively narrow electron scattering profile from the CSM and the damping of the red wing from the scattering from the CDS and the ejecta.

Besides this general agreement with the observations, the model provides several interesting results. One is that the optical depth to electron scattering required to produce the necessary smoothing of the input emission only slowly decreases with time. In Sect. 8.3, we provide possible explanations to this result.

The O I λ11287 line profile is similar to Hα, although the blue extension is somewhat smaller (see Fig. 16). This line may arise as a result of a recombination from O II, collisional excitation, or fluorescence. If recombination dominates, one expects the strength of other high excitation lines to be similar, in particular the O I λ7774 line. This line is, however, very weak if it is detected at all (see Fig. 7), and we therefore exclude this mechanism. Collisional excitation is also unlikely since this would also result in comparable line flux for the O I λ7774 line.

A likely explanation for the excitation is fluorescence by the Lyman-β line. This was discussed by Grandi (1975) for AGNs and has also been invoked for SNe and, in particular, in the case of SN 1987A. This process results in the line O I λ11287 and then O I λ8446, with similar flux values. The latter line is, however, blended with the Ca II triplet as well as a possible strong Fe II line (Fig. 7), but it is consistent with the observations. Therefore these lines are likely to come from the hydrogen-rich ejecta, and not from processed gas. We also note that the calcium lines are likely produced in the same regions as Hα since the Ca II lines can be reproduced using a combination of Hα profiles.

Finally, we comment on the origin of the Hα emission during the first month of SN 2013L. As shown in Fig. 9, the Hα is symmetric until +27 d. The emission from the shock at these early phases is thermalized due to the large optical depth of the CSM. Only at +27 d do the photons from the shock escape without thermalization and the line profile shows asymmetry. Therefore, at very early phases, the Hα emission comes from the unshocked gas, where photons with τ_e ≲ 10 can escape.

7.4. A CSM-interaction model

After ∼20−30 days, we found from the line profile fitting that τ_e ≲ 10. The diffusion time scale is therefore t_diff ≈ τ_eR/c ≈ 7(v/c)t ≈ t/10, so t_diff ≪ t. At times later than a month, it should therefore be a good approximation to ignore diffusion for the light curve. This is also supported by the simulations of Dessart et al. (2015, see their Fig. 3), where the energy input from the shock is compared to the bolometric luminosity.

If radiative and with 100% efficiency (see below), the luminosity from the forward shock is given by

$$\begin{array}{c}\hfill L=2\pi {\rho }_{\mathrm{CSM}}{r}_{\mathrm{sh}}^2{V}_{\mathrm{sh}}^3=\frac{1}{2}\frac{\dot{M}}{u_{\mathrm{wind}}}{V}_{\mathrm{sh}}^3,\end{array}$$(10)

assuming a steady mass-loss rate for the CSM. In Fig. 25, we plotted this luminosity as a function of time, assuming a shock velocity from Eq. (9). For comparison purposes, we also plotted the bolometric luminosity of SN 2010jl from Fransson et al. (2014). In this fit for SN 2010jl, we used the same velocity dependence as for SN 2013L, since a shock velocity could not be derived from the line profiles in that case. The assumed velocity is, however, compatible with that inferred from X-ray observations (Ofek et al. 2014b). For the velocity, we used u_wind = 120 km s⁻¹ for SN 2013L (Sect. 5.1) and u_wind = 105 km s⁻¹ for SN 2010jl (Fransson et al. 2014). As seen in Fig. 25, this simple model gives a good fit to the observations for both SNe up to ∼300 days, after which the observed light curve becomes steeper. We note that the points that depart most from the model are the luminosities derived from using only NIR spectra, which are also the ones affected most by the uncertain extrapolation into the optical range.

Fig. 25. Bolometric luminosity as a function of time for SN 2013L together with SN 2010jl. The dashed lines represent the luminosity predicted for a model with the velocity dependence from Eq. (9) and mass-loss rates Ṁ = 1.7 × 10−2 M⊙ yr−1 and Ṁ = 9 × 10−2 M⊙ yr−1, respectively. For SN 2013L the blue dots refer to luminosities derived from optical spectral fits, while the diamonds refer to NIR fits.

With these assumptions, we derived the mass-loss rates of Ṁ = 1.7 × 10⁻² M_⊙ yr⁻¹ for SN 2013L compared to Ṁ = 9 × 10⁻² M_⊙ yr⁻¹ for SN 2010jl. Because the CSM velocity of these two SNe were very similar, the mass-loss ratio also corresponds to a similar ratio for the densities.

Extrapolating back to the explosion and assuming a shock velocity of V_s(t) = V₀(t₀/t)^α (as in Eq. (9)), the total CSM mass swept up by the shock is

$$\begin{array}{c}\hfill M_{\mathrm{shell}}=\frac{\dot{M}{V}_0{t}_2}{(1-\alpha )u_{\mathrm{wind}}}{\left(\frac{t_0}{t_2}\right)}^{\alpha }·\end{array}$$(11)

Here, t₂ is the time the interaction effectively stops, which we consider to be 350 days when the break in the bolometric light curve is seen. With V₀ = 4000 km s⁻¹, t₀ = 100 days, α = 0.23, and the above values for Ṁ and u_w, we find M_shell = 0.64 M_⊙. The duration of this high mass-loss phase is ∼M_shell/Ṁ ≈ 42 years. A similar duration can be obtained from the drop in luminosity at ∼350 days and the shock velocity at this epoch, t_eruption ≈ t₂V(t₂)/u_w ≈ 350 × 2700/105 = 25 years.

The mass-loss rate found here is larger than the rate found by A17, but it is similar to those found for other SN IIn progenitors. At +33 d, A17 found a mass-loss rate of 0.34−1.75 × 10⁻³M_⊙ yr⁻¹, using $\dot{M}=2u_{\mathrm{wind}}L(\mathrm{H}\alpha )ϵ{(\mathrm{H}\alpha )}^{-1}{V}_{\mathrm{sh}}^{-3}$, where ϵ(Hα) is the fraction of the bolometric luminosity emerging in the Hα line. Their assumed value for this, ϵ(Hα) = 0.1−0.5, is very high compared to what is found here, ϵ(Hα) = 6 × 10⁻³. Also, their velocity, V_sh = 2696 km s⁻¹ (from the electron scattering wings of Hα) is low by a factor of ∼2 compared to what we find. Together, this explains the discrepancy.

A problem with the above mass-loss rate is that the optical depth to electron scattering of the CSM is lower than estimated from the line profiles in Sect. 7.2. Even at +48 days, the optical depth corresponds to only ∼1.7. This is even more accentuated by the fact that the lines are still optically thick at +330 days and even the day +728 observation shows a clear electron scattering wing (Fig. 14). This is close to or after the drop in luminosity, which may be interpreted as the break out of the CSM by the shock and consequently a drop in the CSM column density.

A caveat to the above mass-loss rate and total mass estimates is that multidimensional effects likely affect the radiative efficiency for radiative shocks through thin-shell instabilities, creating a corrugated shock interface (Kee et al. 2014; Steinberg & Metzger 2018). This is found to decrease the ratio between the X-ray luminosity and the total luminosity from the shocks, but it is not clear if this changes the total radiative efficiency. A larger fraction of UV and soft X-rays actually increases the absorbed fraction of the radiation because of the larger photoabsorption cross sections. If the radiative efficiency is lowered by a factor of η, this increases the mass-loss rate and total mass by a factor of η⁻¹. An anisotropic mass-loss rate has a similar effect. For the simple case of a conical outflow with solid angle, Ω, the mass-loss rate inferred from the luminosity increases by a factor of 4π/Ω for a fixed luminosity. Both from observations of LBVs, such as Eta Car (Smith 2006), and from observations of SN 2010jl (Fransson et al. 2014; Patat et al. 2011), for example, there is evidence for asymmetries in the CSM.

A further caveat is that we may have underestimated the bolometric luminosity due to bolometric corrections in the UV and IR, or due to reddening. We have already discussed this in Sect. 6 with regard to the UV, which may be severely underestimated. In addition, some of the X-rays from the shock may escape without being thermalized, depending on the shock velocity.

The cross section for photoelectric absorption can be approximated by σ_abs ≈ σ_T(E/10 keV)⁻³. The effective optical depth for thermalization is ${\tau }_{\mathrm{therm}}\approx \sqrt{{\tau }_{\mathrm{abs}}{\tau }_{\mathrm{e}}}\approx {\tau }_{\mathrm{e}}{(E/10\mathrm{keV})}^{-3/2}$. For τ_e = 5, one finds τ_therm ≈ 1 at ∼30 keV, which is similar to the shock temperature for a shock velocity of 5000 km s⁻¹ (Eq. (1)). The thermalization of the X-rays may, therefore, only be fulfilled for the low energy part of these. As discussed earlier, the shock properties and the radiation are, however, not yet fully understood.

We include these effects in a parameter η, similar to that invoked by Ofek et al. (2014b). A more general form of Eq. (10), for an arbitrary density slope, is then

$$\begin{array}{c}\hfill L(t)=\frac{1}{2}\eta \frac{\dot{M}}{u_{\mathrm{wind}}}{V}_0^3{\left(\frac{t}{t_0}\right)}^{[n(2-s)+6s-15]/(n-s))}.\end{array}$$(12)

We can estimate the velocity, optical depth, and luminosity from the observations. At t₀ = 48 days, we estimated V₀ = 4800 km s⁻¹, τ_e(t₀) = 8 and L(t₀) = 6 × 10⁴² erg s⁻¹.

We assumed that n = 5 and varied the power law index of the density law s, which for this n is only important for the evolution of τ_e. The time of the breakout is not well determined. For a spherical geometry, one may consider this to be the time of the drop in luminosity, ∼350 days. Deviations from this are discussed below. A calculation with these parameters is shown in the left panel of Fig. 26 for s = 0, 1 and 2. The most interesting result is that to reconcile both the luminosity, optical depth, and velocity, a mass-loss rate corresponding to ∼0.15 M_⊙ yr⁻¹ is necessary while the efficiency parameter is low, η ∼ 0.14 for s = 2.

Fig. 26. Models for the luminosity, electron scattering optical depth of the CSM (solid lines) and CDS (dashed lines), and velocity based on the similarity solution for different values of the density power law index s. Left panels: break out time tb = 350 days; right panels: tb = 700 days. The mass-loss rate and efficiency parameters are given in the legends for the different s.

As shown in the second panel, the sharp drop in density at 350 days, however, results in a very low value of τ_e, which is in contradiction to the optical depth at 330 days (Sect. 7.2). What we believe to be the most likely way of reconciling both a value of τ_e ≳ 1 and a drop in luminosity would be if the CSM is asymmetric, as discussed above.

An asymmetric CSM would result in different times for the shock break-out and a steepening of the light curve, as shown in van Marle et al. (2010) and Vlasis et al. (2016). The optical depth of the CSM in the directions of large column density may, however, still be large corresponding to a larger value of R_out, that is, t_b and τ_e in this direction (Eqs. (4)−(12)). In reality, an asymmetric CSM would also affect the light curve and velocity. To model this, a multidimensional hydrodynamic calculation would be necessary, as discussed in van Marle et al. (2010) and Vlasis et al. (2016). To illustrate the effect on τ_e, we show one calculation using t_b = 700 days. The longer t_b results in a less steep decrease of τ_e, as expected, and a lower mass-loss rate (because R_out increases). For s = 0, the value of τ_e is nearly constant. To get the same luminosity, the efficiency has to increase (see Fig. 26).

As seen from (Eqs. (4)–(12)), the value of R_out, that is, t_b, only affects τ_e. In reality, an asymmetric CSM would also affect the light curve and velocity. To model this, a multidimensional hydrodynamic calculation would be necessary, as discussed in van Marle et al. (2010) and Vlasis et al. (2016). To illustrate the effect on τ_e, we show one calculation using t_b = 700 days. The longer t_b results in a less steep decrease of τ_e, as expected, and a lower mass-loss rate (because R_out increases). For s = 0, the value of τ_e is nearly constant. To get the same luminosity, the efficiency has to increase (see Fig. 26).

8. Discussion

8.1. The progenitor system

A first indication of the nature of the progenitor system of SN 2013L can be obtained from the observed narrow Hα component, which provides a direct measure of the unshocked CSM wind velocity. The origin of this CSM is likely to be the stellar winds of the progenitor that enriched the environment over its evolutionary lifetime. A wind velocity around 120–240 km s⁻¹ is consistent with a luminous blue variable wind (see e.g., Kiewe et al. 2012; Taddia et al. 2013; Smith 2014), but also with a wind from a yellow hypergiant (YHG; Smith 2014; Smith et al. 2015), as is also noted by A17. Typical red supergiant winds have lower velocities (∼10 km s⁻¹; Smith et al. 2009), whereas typical Wolf-Rayet wind velocities are higher (≳1000 km s⁻¹; Pastorello et al. 2007).

The mass-loss rate is at least about 1.7 × 10⁻²M_⊙ yr⁻¹ and from our discussion in the previous section, it is more likely to be ∼0.15 M_⊙ yr⁻¹. For these respective mass-loss rates, this is consistent with an LBV progenitor that is capable of expelling a large amount of mass, from 0.4−0.7 M_⊙ up to 3.8−6.3 M_⊙ in 25−42 yr in this case. The most well-studied case of this is η Carinae, which expelled about 10 M_⊙ in 30 yr during the great eruption that occurred in the 19th century (Smith 2013). Woosley (2017) discusses this case in the context of pulsational pair instability SNe. We note that in the case of SN 2013L, there is no direct indication of advanced nuclear burning, which is as expected for this type of event. A YHG progenitor can be considered if the mass-loss rate is closer to 10⁻³M_⊙ yr⁻¹ rather than 10⁻²M_⊙ yr⁻¹.

Our spectral model does not require strong asymmetry of the CSM to reproduce the Hα profile, but rather a combination of occultation and electron scattering. However, we cannot exclude the possibility of asymmetry in the CSM, as in the case of η Carinae and its complex CSM characterized by polar lobes (see, e.g., Smith 2013). In fact, as we discuss in Sect. 7.4, there are indications in favor of this from the light curve and line profiles.

The large CSM mass indicated by our analysis is also consistent with the strong deceleration of the Hα emitting CDS, with V_shell ∝ t^−0.23 (Fig. 24, Eq. (9)), showing that a large fraction of the kinetic energy of the ejected material is radiated. Such, or even more extreme scenarios have been invoked by Dessart et al. (2016) and Chugai (2016) for SN 1994W and SN 2011ht. In these simulations, the ejecta mass is much lower than the CSM mass, resulting in, in some cases, complete deceleration of the fast ejecta and an initially luminous SN, dropping fast at ∼100 days. In the case of SN 2013L, the velocity is still ≳1700 km s⁻¹ at ∼1500 days (Fig. 24) and the drop in luminosity is slower than in SN 1994W. However, the need for a very massive CSM is probably required, although the ejecta mass may be higher than for SN 1994W.

8.2. Dust in SN 2013L

In SN 2013L, there are clear signs, at least at late epochs, of luminous NIR and MIR emission, which is commonly associated with dust emission. The dust could already be present in the unshocked CSM and from there it could easily drive a light echo where the SN optical and UV emission is redistributed to longer wavelengths, but it could also be associated with newly formed dust in the SN ejecta (e.g., Stritzinger et al. 2012; Gall et al. 2014; Chugai 2018). In the case of SN 2010jl, an origin in the inner ejecta has been argued against (Fransson et al. 2014; Sarangi et al. 2018). Here we attempt to determine where the dust is located in SN 2013L and characterize the nature of the NIR emission.

The NIR and MIR emission was modeled by fitting a single BB function to the NIR and MIR part of the SED that corresponds to the warm component discussed in Sect. 6. We find that between +132 d and +861 d, the luminosity (∼10⁴² erg s⁻¹), the temperature (1300 K), and the radius (10¹⁶ cm) of this warm component are rather constant.

In Fig. 27, we plotted the radius of the warm BB component already shown in Fig. 20. If produced by spherically-distributed dust with a covering factor (Ω_CF) of 4π, then the BB radius is approximately the dust radius. If the dust is distributed with a smaller covering factor (see the cases of SNe 2005ip and 2006jd, Stritzinger et al. 2012), the dust radius is larger than the BB radius. If the dust is located at a certain radius from the SN, it is necessary check if the forward shock could reach and destroy the dust. We adopted the velocity of the shell from Fig. 23 as a measure of the velocity of the forward shock, which evolves as a power law of index −0.23. The radius reached by the forward shock (blue lines) intercepts that of the dust already at ∼+300 d if we assume Ω_CF = 4π (red lines and squares). For a smaller covering factor, Ω_CF = 4π/5, the dust is not reached by the forward shock, and it can survive and lead to echo emission. We also checked whether or not the high SN luminosity could vaporize the dust using the technique adopted in Stritzinger et al. (2012). For a covering factor Ω_CF = 4π/5, the dust should consist of grains larger than ∼0.01 μm in order to survive. For a Ω_CF = 4π, the dust would be vaporized by the high SN luminosity, regardless of its grain size. We find that it is entirely plausible that the dust in SN 2013L is pre-existing and it covers only a fraction of the 4π surface.

Fig. 27. Radius of the warm BB component (NIR emission) assuming it originated from dust distributed spherically with full covering factor (red) and with a smaller covering factor (4π/5, magenta). The lines are from the SEDs fit and the symbols are from the spectral fit with BBs. The radius of the forward shock as obtained from the spectral modeling of Hα (see Fig. 23) is shown in blue. The velocity of the reverse shock from which we derived its radius (shown in green) is taken to be 1000 km s−1 slower than that of the forward shock. The evaporation radius for the dust, which depends on the dust grain size (see legend), is shown in black. If the NIR emission is due to pre-existing dust, this must be distributed with a covering factor of ≲4π/5 otherwise the ejecta would have already reached and destroyed it around 300 days. Furthermore, the grain size of the dust must be <0.1 μm otherwise the large luminosity of the SN would have vaporized it (this is assuming a covering factor of 4π/5).

Given the uncertainty on the dust radius, there might also be the possibility that the dust is formed, with Ω_CF = 4π, within the cold dense shell between the forward shock and the reverse shock, which we have assumed to move 1000 km s⁻¹ slower than the forward shock, see green lines in Fig. 27. However, we do not see an increasing depression on the red side of the emission lines, such as Hα, which would be a sign of newly formed dust. In Stritzinger et al. (2012), SN 2005ip and SN 2006jd show this progressive depression on the red side of Hα at relatively early epochs. There is also the possibility that the dust did not form close to the line forming regions, and therefore we did not observe the depression. Furthermore, the high dust temperature at late epochs might point to dust formation in the ejecta at these late phases (see Gall et al. 2014; Chugai 2018, for the case of SN 2010jl).

8.3. The shock structure

In Fig. 25, we show the bolometric light curves of SN 2010jl and SN 2013L. Although the luminosity of SN 2010jl was a factor of ∼5 higher, the general shape is similar with a break ∼300−400 days and a similar slope of the light curve, both before and after the break. While the paper by Dessart et al. (2015) aimed to model SN 2010jl to a large extent, many of the results are therefore also relevant to SN 2013L.

From the modeling of the line profiles, we find several interesting results, which give information about the structure of the shock and dense CSM. A unique feature of SN 2013L is that we have direct information about the shock velocity from the skewed line profile. One may then ask what the difference is between SN 2013L and other Type IIn SNe, which show more symmetric line profiles. As Fig. 28 shows, this is mainly a result of a lower optical depth to electron scattering of the dense unshocked CSM from the shock. In this figure, we show the line profile for several optical depths, including the one used for the fit at day 48 in Fig. 24τ_e = 12.3. As in these calculations, we used a maximum electron temperature T_e = 1.5 × 10⁴ K. One can obtain the same line profiles for other temperature profiles, as seen from scaling ${\tau }_{\mathrm{e}}\propto T_{\mathrm{e}}^{-1/2}$ for a constant temperature medium.

Fig. 28. Line profiles for different optical depth of the CSM. The black with τe = 7.9 corresponds to the best fit model at +48 days. The velocity of the emitting shell behind the shock is in all cases Vem in = 5390 km s−1 and Vem out = 6160 km s−1. One notes the gradual transition from highly skewed profiles to symmetric as the optical depth increases. In all cases, an electron temperature of 15 000 K is assumed.

As we see from Fig. 28, once τ_e ≳ 10, the shoulder in the line profile disappears and the line becomes increasingly symmetric as the optical depth increases further. The most likely explanation for the usually symmetric lines is, therefore, that for most other bright Type IIn SNe, the shock has been at such a large optical depth to electron scattering that the velocity information has been lost.

Another special feature of SN 2013L is that we can identify two different components, one central, symmetric and one broad, skewed profile. We interpret these to come from the preshock and postshock regions, respectively. As is seen in Fig. 28, for low optical depths (e.g., the τ_e = 3.75 profile), these components may give rise to a double peak structure, one narrow component from the preshock gas, and a broader component from the postshock. At the other extreme, we note the almost perfectly symmetric line for τ_e = 30.8. This is in spite of the high shock velocity and also the occultation by the photosphere, which is assumed to coincide with the inner boundary of the emitting region, here at V_ph = V_em in = 5390 km s⁻¹. At lower optical depths, the occultation explains the depression on the red side of the line.

A somewhat surprising result that was found is that the optical depth to electron scattering required to produce the necessary smoothing of the input emission is only weakly dependent on time, τ_e = 12.3 at +48 d, compared to τ_e = 7.2 at +330 d. The difference is hardly significant and depends on the assumed temperature of the gas, for example. Naively, one could expect this to decrease by a large factor as the column density of the CSM in front of the shock decreases. As discussed in Fransson et al. (2014, see their Eq. (17)), if ionization and not density bounded, the optical depth of the region ionized by the shock radiation is for a an ρ ∝ r⁻² CSM mainly sensitive to the shock velocity and not to the density of the CSM. The lower density results in a lower shock luminosity for a radiative shock ($\propto n_{\text{CSM}}{V}_{\text{s}}^3$), but also in a lower total recombination rate ($\propto n_{\text{CSM}}^2\Delta r$). The optical depth, τ_e ∝ n_CSMΔr, of the ionized region is therefore sensitive to the shock velocity, ${\tau }_{\text{e}}\propto V_{\text{s}}^3$, but independent of density. Because V_s is only slowly decreasing with time, the electron scattering depth is expected to decrease slowly.

For SN 2010jl, Dessart et al. (2015) find τ_e ≈ 5−10 between 30 and 100 days, decreasing to unity at ∼350 days. The fact that this model shows a decrease in the optical depth probably means that the ionized region is density bounded.

The ${\tau }_{\text{e}}\propto V_{\text{s}}^3$ scaling may also explain why Type IIn SNe in many cases have an electron scattering depth so large that we do not see any velocity shift. As discussed in Sect. 7.4, an anisotropic, dense CSM may resolve the problem with the drop in the light curve, but τ_e that is still high.

9. Conclusion

To conclude our analysis, we schematize the most likely CSM interaction scenario for SN 2013L in Fig. 29. Here the ejecta of a LBV (or YHG) star interacts with its CSM and thereby produces both a forward and a reverse shock. The spectral model provides a good fit of the Hα profile without invoking asymmetry of the CSM, which we represent as spherically symmetric, which is likewise for the ejecta. The high luminosity of the SN destroys all the dust within the evaporation radius. Outside the evaporation radius, pre-existing dust reprocesses the UV and optical light of the SN into the NIR and MIR excess that we observe. The dust is distributed with a small covering factor, which might hint toward some asymmetry in the outer CSM. The Hα broad component (red line in the Hα inset) is produced in the cold-dense shell between forward and reverse shock (dashed area) and the red side of this emission is occulted by the ejecta. The broad, boxy emission of the shell (blue line in the Hα inset) is broadened by Thomson scattering due to CSM located outside the shock. The central component (cyan in the Hα inset) is produced by narrow emission in the CSM directly outside the forward shock which is again broadened by Thomson scattering in the CSM. The excited, unshocked, outer CSM produces the narrow Hα component, which reveals a slow wind (120−240 km s⁻¹) that is compatible with a LBV (or a YHG) star. Our mass-loss-rate estimates favor a LBV progenitor star. We find strong arguments from the light curve and electron scattering wings for an anisotropic CSM.

Fig. 29. Schematic view of SN 2013L interacting with its CSM, as inferred from our analysis.

Acknowledgments

We thank the referee for their comments which helped to improve the paper, especially concerning the modeling. We thank Jennifer Andrews for having shared with us the late-time spectra of SN 2013L. F.T. and M.D.S. acknowledge support for a project grant provided by the Independent Research Fund Denmark (IRFD). F.T. and M.D.S. also acknowledge support from the VILLUM FONDEN under experiment grant 28021. F.T. and J.S. gratefully acknowledge the support from the Knut and Alice Wallenberg Foundation. C.F. acknowledges support for the Swedish Research Council. M.D.S. acknowledges support from a research grant (13261) from VILLUM FONDEN. The Oskar Klein Centre is funded by the Swedish Research Council. The CSP-II has been supported by the National Science Foundation under grants AST1008343, AST1613426, AST1613455, and AST1613472, and also in part by a grant from the Danish Agency for Science and Technology and Innovation through a Sapere Aude Level 2 grant (PI M.D.S.). The Swift Optical/Ultraviolet Supernova Archive (SOUSA) is supported by NASA’s Astrophysics Data Analysis Program through grant NNX13AF35G. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Additional tables

All Tables

All Figures

	Fig. 1. Finding chart for SN 2013L, from a combination of all of the BVgri images obtained with the Swope telescope between JD 2456316.7 and JD 2456448.6.
In the text

Fig. 2. CSP-II uBgVriYJHKs-band and Swiftuvw2, uvm2, uvw1-band light curves of SN 2013L (filled symbols). Light curves were shifted for clarity by the amounts shown in the legend. The photometry is also included from A17 (open symbols). Solid lines correspond to tension spline fits to the CSP-II photometry. Epochs of spectral observations are indicated with vertical segments (NIR spectra at the top and optical spectra at the bottom; in black our new spectra, in red those already published by A17). The nondetection limit in an unfiltered/r band is indicated by a red triangle, while unfiltered/r band discovery and confirmation magnitudes are shown by red stars. The top-right inset shows a zoomed in view of the light curves over the first 150 days of evolution.

In the text

Fig. 3. Visual-wavelength spectroscopic time-series of SN 2013L. Indicated next to each spectrum is the phase, i.e., days from discovery. The spectra that were already published in Andrews et al. (2017) and publicly available are marked by a star next to their phases, the four late-time spectra from A17 are marked by two stars. Each spectrum was normalized by its median and shifted by a constant for better visualization.

In the text

Fig. 4. NIR spectroscopic time-series of SN 2013L obtained with FIRE and X-shooter. The phase of each spectrum relative to discovery is indicated next to each spectrum. The spectra that were already published in Andrews et al. (2017) are marked by a star next to their phases. Each spectrum was normalized by its median and shifted by a constant for a better visualization. We marked the telluric bands at the top.

In the text

	Fig. 5. Absolute magnitude light curves of SN 2013L covering UV to NIR wavelengths as well as those of the long-lasting Type IIn SN 1988Z (Turatto et al. 1993), SN 1996al (Benetti et al. 2016), SNe 2005ip and 2006jd (Stritzinger et al. 2012), SN 2010jl (Fransson et al. 2014), iPTF13z (Nyholm et al. 2017), and KISS15s (Kokubo et al. 2019).
In the text

	Fig. 6. B − V, g − r, and i − H color evolution of SN 2013L as compared to that of other long lasting SNe IIn (SN 1988Z from Turatto et al. 1993, SN 1996al from Benetti et al. 2016, SNe 2005ip and 2006jd from Stritzinger et al. 2012, SN 2010jl from Fransson et al. 2014, iPTF13z from Nyholm et al. 2017, KISS15s from Kokubo et al. 2019).
In the text

Fig. 7. Visual (top) and NIR (bottom) spectral line identifications for SN 2013L. It is important to note that the spectra in the top panel were smoothed (dark lines) for presentation purposes. In the optical, we plotted the expected positions for the Fe II lines (red) for a model with Fe II Lyα and Lyβ pumping, from Sigut & Pradhan (2003). These lines can explain the blue excess in the spectrum at late epochs.

In the text

Fig. 8. Left panel: Hα, Hβ, and Pβ line profiles compared to one another and plotted in velocity space. Each line has been corrected for extinction and continuum subtracted. Hβ and Pβ lines were scaled by the constants indicated on the left to match the intensity of Hα. The spectral phases in days since discovery are reported on the right. The shape of the line profiles is remarkably similar at all epochs. Top-right panel: Hα luminosity and continuum luminosity in the Hα region as a function of time for SN 2013L and SN 2010jl (Fransson et al. 2014). Bottom-right panel: Hα/Hβ and Hβ/Hγ flux ratio as a function of time for SN 2013L and SN 2010jl.

In the text

Fig. 9. Left panel: continuum-subtracted and extinction corrected Hα profiles (black) fitted by a sum of a Lorentzian (blue) and a Gaussian (green) function to reproduce the broad components. The total fit is shown in red. The spectral phases are shown on the right next to each spectrum. The wavelength range around the narrow component was excluded from the fit. Right panel: velocities and luminosities of the Hα components. The centroid and the FWHM of the Gaussian and of the Lorentzian are shown in the top panels, their luminosities and the total luminosity, as well as the velocities of the narrow component are shown in the bottom panels. Here for the narrow component, we considered the velocity of the absorption minimum assuming zero velocity corresponding to the emission peak, such as in A17. We report our results as well as those of A17 to show the close agreement.

In the text

Fig. 10. Left panel: normalized spectra in velocity space around the [O I] λ6300 wavelength regions. The original spectrum in black and the smoothed spectrum in red are shown to better guide the eye to follow the line profile evolution. The spectral phases are reported in red next to each spectrum. The emergence of the oxygen line at 6300 Å at late epochs is evident. We note that the line at −3000 km s−1 is Fe II. Right panel: continuum subtracted and extinction corrected spectra in velocity space around the wavelength region of O I λ11827. This line appears at late epochs and it is well fit by a single Gaussian profile of ∼4000 km s−1 FWHM (shown in red). The phases of each spectrum are reported next to them.

In the text

Fig. 11. Continuum-subtracted and extinction corrected Ca II triplet profiles (black) fitted by a sum (blue) of three Hα best-fit profiles from Fig. 9 to reproduce the complex line shape. We followed the approach by A17 and confirmed that this is also a good fit for the Ca II triplet at +302 d and +330 d. The spectral phases are shown on the right next to each spectrum. An excess blueward of the best line fit is compatible with the emerging of O I λ8446, but there is also Fe II that can contribute to the flux in that range.

In the text

Fig. 12. Narrow component of Hα plotted in velocity space. The spectral phases in days since discovery are reported on the right. The zero velocity was chosen to correspond to 6562.8 Å in the rest frame. We plotted (in red) the best fit to the line profile, which was achieved by fitting two Gaussian components, one for the absorption and one for the emission. The emission dominates at early epochs, while absorption becomes stronger at later epochs. We marked the BVZI with a dashed segmented line that ranges from −120 to −240 km s−1.

In the text

	Fig. 13. Evolution of the continuum subtracted Hα line profile in velocity space. The flux was normalized to the flux of the flat part of the profile between the shoulder and the central component. The spectra from +482 d to +1133 d are from A17.
In the text

Fig. 14. Continuum-subtracted and extinction corrected Hα profiles (black), this time fitted by a sum of a skewed, broad Gaussian (blue) and a narrower Gaussian (green) function to reproduce the broad components. The total fit is shown in red. The spectral phases are shown on the right next to each spectrum. The wavelength range around the narrow component was excluded from the fit. Right panel: velocities and luminosities of the Hα components. The centroid and the FWHM and HWHM of the skewed broad Gaussian (blue) and of the central Gaussian (green) are shown in the top panels, their luminosities and the total luminosity, as well as the velocity of the narrow component are shown in the bottom panels. Here, for the narrow component, we consider the blue-velocity-at-zero-intensity (BVZI) assuming zero velocity at the Hα rest wavelength to obtain a better estimate of the CSM velocity.

In the text

	Fig. 15. Flux as a function of time for three bright hydrogen lines in the NIR spectra. Their evolution is rather similar.
In the text

	Fig. 16. Continuum-subtracted Hα and O I λ11287 at +300 to +330 d. The NIR O I line, which was scaled to match the peak of Hα, shows a broad profile that is almost as broad as Hα.
In the text

Fig. 17. Inset: optical spectra after absolute flux calibration and extinction correction, fitted by a single BB function. The parts of the spectrum in light gray were excluded from the fit because of the presence of strong emission lines. Main panel: NIR spectra after absolute flux calibration (using J band) and extinction correction, fitted either by a single BB (at early epochs), or by a sum (magenta) of two BB components, one warm (red) and one hot (blue). The warm (red) component dominates at late epochs. The flux of the last spectrum was multiplied by five for a better visualization.

In the text

Fig. 18. Example of SEDs of SN 2013L with the specific epochs corresponding to times when NIR photometry is available. Each flux point was obtained by converting the magnitudes into fluxes after having applied the appropriate extinction correction. The color and symbol code for each band is the same as in Fig. 2. A single BB function was fit to each SED until +68 d (black line), while a sum of two BB functions was used to fit the +335 d SED in order to reproduce the NIR excess (red line) and the optical BB shape (blue line). The r band is not included in the fit due to the prevalent Hα emission. The BB temperature is reported in each subpanel as well as the phase of the SED.

In the text

Fig. 19. Late phase Ks-band and Spitzer photometry converted into monochromatic fluxes plotted together with our last epoch NIR spectrum. The SED exhibits a strong NIR and MIR excess that dominates the SN emission two years after discovery. The NIR spectrum was scaled to match the Ks-band flux. We fit a sum of two BB functions (magenta dashed line) to reproduce the NIR spectrum and the MIR fluxes from the first Spitzer epoch. The warm BB component (red dashed line) dominating the NIR and MIR has a temperature of 1022 K and the hot component (blue dashed line) has a temperature of 4513 K.

In the text

Fig. 20. Top panel: bolometric luminosity of SN 2013L. From the BB fit to the photometry (solid lines), we derived the luminosity of the hot BB component (solid blue line, dominating at early epochs) and of the warm BB component (red solid line, dominating at late epochs and emerging at +132 d). Their sum (dashed magenta line) and their sum plus the Hα emission neglected from the BB fit (black solid line) are also shown, the last one represents the total luminosity. We also obtained the corresponding luminosities from the BB fit to the optical (blue diamonds) and NIR spectra (blue dots and red squares). Central panel: temperature evolution of SN 2013L for the two different BB components fitting the SEDs of SN 2013L (symbols as in the top panel). In addition, the dashed lines mark the uncertainties on the temperature from the BB fit to the photometry. Bottom panel: radius evolution for the two different BB components fitting the SEDs of SN 2013L, symbols are the same as those in the top and central panels.

In the text

	Fig. 21. Fit of the radius evolution used to derive the explosion epoch. A linear fit from the first 10 days (red) gives an explosion epoch of −15 d before explosion. A similar explosion epoch was obtained by fitting the first +30 d with a second-order polynomial (magenta).
In the text

Fig. 22. Continuum subtracted Hα line (shown in black) from days 48 (left) and 330 (right) together with the electron scattering profiles from the model with Hα emission only from the ejecta and shocked region (red). The dashed magenta lines show the box-like input emission profiles from a radially thin shell with velocity 4800 km s−1 and 2700 km s−1, respectively. The cyan lines show the residual between the model and the observed profile. Electron scattering and occultation by the photosphere result in the smooth, blue-shifted profile. We note the central excess from the models, indicating a different component.

In the text

	Fig. 23. Red squares: velocity evolution of the Hα-line blue shoulder as function of time. Blue squares: velocity of the shell from the electron scattering model for selected dates. The dashed lines show a least square fit to these.
In the text

Fig. 24. Same as Fig. 22, but now also including Hα emission from the unshocked medium preionized by the shock, which can be seen as the narrow magenta line on top of the broad from the shocked gas (normalized to one at maximum). The magenta, dashed profile shows the line profile without electron scattering or continuum thermalization. We note that the residual close to the line center is now much lower at both epochs. The steep profile on the red side is a result of the comparatively narrow electron scattering profile from the CSM and the damping of the red wing from the scattering from the CDS and the ejecta.

In the text

Fig. 25. Bolometric luminosity as a function of time for SN 2013L together with SN 2010jl. The dashed lines represent the luminosity predicted for a model with the velocity dependence from Eq. (9) and mass-loss rates Ṁ = 1.7 × 10−2 M⊙ yr−1 and Ṁ = 9 × 10−2 M⊙ yr−1, respectively. For SN 2013L the blue dots refer to luminosities derived from optical spectral fits, while the diamonds refer to NIR fits.

In the text

	Fig. 26. Models for the luminosity, electron scattering optical depth of the CSM (solid lines) and CDS (dashed lines), and velocity based on the similarity solution for different values of the density power law index s. Left panels: break out time tb = 350 days; right panels: tb = 700 days. The mass-loss rate and efficiency parameters are given in the legends for the different s.
In the text

Fig. 27. Radius of the warm BB component (NIR emission) assuming it originated from dust distributed spherically with full covering factor (red) and with a smaller covering factor (4π/5, magenta). The lines are from the SEDs fit and the symbols are from the spectral fit with BBs. The radius of the forward shock as obtained from the spectral modeling of Hα (see Fig. 23) is shown in blue. The velocity of the reverse shock from which we derived its radius (shown in green) is taken to be 1000 km s−1 slower than that of the forward shock. The evaporation radius for the dust, which depends on the dust grain size (see legend), is shown in black. If the NIR emission is due to pre-existing dust, this must be distributed with a covering factor of ≲4π/5 otherwise the ejecta would have already reached and destroyed it around 300 days. Furthermore, the grain size of the dust must be <0.1 μm otherwise the large luminosity of the SN would have vaporized it (this is assuming a covering factor of 4π/5).

In the text

Fig. 28. Line profiles for different optical depth of the CSM. The black with τe = 7.9 corresponds to the best fit model at +48 days. The velocity of the emitting shell behind the shock is in all cases Vem in = 5390 km s−1 and Vem out = 6160 km s−1. One notes the gradual transition from highly skewed profiles to symmetric as the optical depth increases. In all cases, an electron temperature of 15 000 K is assumed.

In the text

	Fig. 29. Schematic view of SN 2013L interacting with its CSM, as inferred from our analysis.
In the text