Supernova 2014C: Ongoing Interaction with Extended Circumstellar Material with Silicate Dust

2.1.1. Near-IR 1–2.5 μm Photometry

We obtained five epochs of photometry of SN 2014C in the 1–2.5 μm JHKs bands using the Wide-field InfraRed Camera (WIRC; Wilson et al. 2003, with a detector upgrade described in Tinyanont et al. 2019) on the 200 inch Hale Telescope at Palomar Observatory (P200 hereafter). Not every epoch had all three bands, and some epochs had different bands taken a few days apart. NGC 7331 was small enough to fit into the field of view of WIRC, so the data were taken with dithering patterns that send the galaxy around the field of view to measure and subtract the sky background. Science images were dark subtracted and flattened using a flat-field image obtained from median combining dithered sky images with sources masked out. The photometric zero point was determined using ∼30 stars in the field of view with magnitudes from the Two Micron All-Sky Survey (2MASS; Milligan et al. 1996; Skrutskie et al. 2006). The data reduction was performed using an automated IR imaging pipeline described in De et al. (2019).

Figure 1 shows a false-color image using WIRC images from 2018 June 23 where red, green, and blue correspond to the Ks, H, and J bands, respectively. Five images in the right of the figure are from the three near-IR bands and the two Spitzer/IRAC bands to visually demonstrate the red color of SN 2014C at this epoch. As shown in the image, SN 2014C was located in a spiral arm of NGC 7331 with a bright and spatially varying background, making host subtraction difficult. We did not have a pre-explosion image of the host galaxy at a comparable resolution to conduct image subtraction. As a result, we measured the near-IR photometry as follows. First, for each band in each epoch, we created a library of point-spread function (PSF) using stars in the vicinity of the SN that were not on top of the host galaxy. The PSF construction was done using the EPSFBuilder module in the photutils package (Bradley et al. 2019). We used the DAOStarFinder module to find the (subpixel) location of the SN. The PSF photometry was obtained by iteratively subtracting the PSF with different total flux to minimize the summed squared gradient of the image. This quantity was used as the metric to indicate the "point sourceness" of the image as the slowly varying background of the galaxy had a lower gradient than the point source. We compared the residual image of the same band from different epochs to ensure that the SN light had been consistently subtracted for all epochs. The PSF subtraction routine described above was custom-built in Python. Along with these new photometry, we also plotted near-IR photometry published in T16 in Figure 3. All photometry discussed here, and in the following subsections, is listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  IR Photometry of SN 2014C

Date	Epoch	FJ	${\sigma }_{{F}_{J}}$	FH	${\sigma }_{{F}_{H}}$	${F}_{{K}_{s}}$	${\sigma }_{{F}_{{K}_{s}}}$	F[3.6]	${\sigma }_{{F}_{[3.6]}}$	F[4.5]	${\sigma }_{{F}_{[4.5]}}$	Instrument
(UT)	(day)	(mJy)
2014-02-19	37.8	⋯	⋯	⋯	⋯	⋯	⋯	4.55	0.08	6.24	0.09	Spitzer/IRAC
2014-09-05	235.8	⋯	⋯	⋯	⋯	⋯	⋯	2.46	0.10	3.17	0.06	Spitzer/IRAC
2014-10-07	267.4	⋯	⋯	⋯	⋯	⋯	⋯	2.92	0.09	3.56	0.06	Spitzer/IRAC
2015-09-19	614.8	⋯	⋯	⋯	⋯	⋯	⋯	3.57	0.09	5.15	0.07	Spitzer/IRAC
2015-09-25	620.8	⋯	⋯	⋯	⋯	⋯	⋯	3.31	0.09	4.54	0.07	Spitzer/IRAC
2015-10-09	634.8	⋯	⋯	⋯	⋯	⋯	⋯	3.39	0.09	5.25	0.06	Spitzer/IRAC
2016-02-17	765.4	⋯	⋯	⋯	⋯	⋯	⋯	2.57	0.01	4.22	0.02	Spitzer/IRAC
2016-02-24	772.7	⋯	⋯	⋯	⋯	⋯	⋯	2.47	0.01	4.16	0.02	Spitzer/IRAC
2016-03-17	794.3	⋯	⋯	⋯	⋯	⋯	⋯	2.43	0.01	4.15	0.02	Spitzer/IRAC
2016-09-18	979.2	⋯	⋯	⋯	⋯	⋯	⋯	2.05	0.02	3.64	0.02	Spitzer/IRAC
2016-10-12	1003.5	⋯	⋯	⋯	⋯	⋯	⋯	2.04	0.02	3.61	0.02	Spitzer/IRAC
2017-02-20	1134.4	⋯	⋯	⋯	⋯	⋯	⋯	1.91	0.02	3.45	0.02	Spitzer/IRAC
2017-07-11	1275.5	⋯	⋯	⋯	⋯	0.37	0.08	⋯	⋯	⋯	⋯	P200/WIRC
2017-07-13	1277.5	⋯	⋯	0.20	0.02	⋯	⋯	⋯	⋯	⋯	⋯	P200/WIRC
2017-09-27	1353.7	⋯	⋯	⋯	⋯	⋯	⋯	1.83	0.01	3.12	0.01	Spitzer/IRAC
2017-10-03	1359.4	0.11	0.02	⋯	⋯	0.44	0.11	⋯	⋯	⋯	⋯	P200/WIRC
2017-11-25	1412.2	0.16	0.02	0.16	0.03	⋯	⋯	⋯	⋯	⋯	⋯	P200/WIRC
2018-02-28	1507.5	⋯	⋯	⋯	⋯	⋯	⋯	1.65	0.02	3.08	0.01	Spitzer/IRAC
2018-06-23	1622.5	0.08	0.02	0.14	0.02	0.39	0.07	⋯	⋯	⋯	⋯	P200/WIRC
2018-06-24	1623.0	⋯	⋯	⋯	⋯	⋯	⋯	1.84	0.14	3.38	0.43	Gemini/NIRIa
2018-07-24	1653.5	0.08	0.02	⋯	⋯	0.40	0.03	⋯	⋯	⋯	⋯	P200/WIRC
2018-08-02	1662.3	⋯	⋯	0.14	0.03	⋯	⋯	⋯	⋯	⋯	⋯	P200/WIRC
2019-03-09	1881.8	⋯	⋯	⋯	⋯	⋯	⋯	1.54	0.02	2.58	0.02	Spitzer/IRAC
2019-04-17	1920.8	⋯	⋯	⋯	⋯	⋯	⋯	1.44	0.02	2.64	0.02	Spitzer/IRAC
2018-06-28	1627.5							F[9.7] = 16		${\sigma }_{{F}_{[9.7]}}=3$		Subaru/COMICS

Note.

^aGemini/NIRI's filters are L' and M ', which are different from Spitzer/IRAC's filters.

Download table as:  ASCIITypeset image

2.1.2. Spitzer/IRAC Photometry

2.1.3. Ground-based 3–5 μm Photometry

SN 2014C's brightness at 3–5 μm permitted L' (3.43–4.13 μm) and M' (4.55–4.79 μm) band photometry from the ground. We obtained those observations with the Near InfraRed Imager and spectrograph (NIRI; Hodapp et al. 2003) on the Gemini North Telescope as part of a fast turnaround program (GN-2018A-FT-108; PI Tinyanont). In these thermal IR bands, the observations were background limited and we used the f/32 camera, deep-well detector bias, and fast read-out mode to maximize the observing efficiency. Raw images were dark subtracted and field flattened with our Python script. The source was observed with dithering patterns that sent the source in and out of the field of view, allowing us to subtract the bright thermal background from the sky and the telescope. The total on-source integration time for the L' and M' bands were 135.15 and 1604 s, resulting in signal-to-noise ratios (S/Ns) of 13 and 8, respectively. We observed a photometric standard HD 203856 (L' = 6.871, M' = 6.840; Leggett et al. 2003) after the target for flux calibration. The observing strategy for the standard star was similar to that of the SN, except that we only dither the star within the field of view. The images in the L' and M' bands are shown in Figure 2 (left and center). We obtained ground-based 3–5 μm photometry in addition to the Spitzer/IRAC to make sure that we have the 3–5 μm coverage concurrent with the N-band imaging described next.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1.4. N-band Imaging with Subaru/Cooled Mid-infrared Camera and Spectrometer (COMICS)

We obtained near-IR 1–2.5 μm spectra of SN 2014C in seven epochs spanning from 1 to 5 yr post-maximum. Table 2 summarizes all the spectra taken. We used the twin medium-resolution, long-slit, echellette spectrographs: TripleSpec on P200 (Herter et al. 2008) and the Near-Infrared Echellette Spectrometer (NIRES) on the Keck telescope,¹⁰ both of which provided simultaneous 1–2.5 μm spectra of a single source. The only difference between TripleSpec and NIRES was the fixed slit widths of 1'' versus 055. The data were taken by either dithering the source along the slit or alternating between the source and the sky to enable sky subtraction. Data from both TripleSpec and NIRES were reduced using a version of Spextool (Cushing et al. 2004), specifically for TripleSpec and NIRES. The software applied field flattening, retrieved a wavelength solution from sky lines present in science observations, and subtracted each exposure pair to remove most of the sky emission. It located the trace of the object in the reduced 2D images, then fit a low-order polynomial across the trace to estimate host background. Finally, the spectral trace was reduced using an optimal extraction algorithm to enhance the S/N while rejecting bad pixels and cosmic-ray hits (Horne 1986). The telluric and flux calibrations were performed using xtellcor with the standard observations of an A0V star either before or after the SN observation (Vacca et al. 2003).

In addition, we also used the Multi-Object Spectrometer for Infra-Red Exploration (MOSFIRE; McLean et al. 2012) on the Keck telescope. MOSFIRE can obtain spectra of up to 46 sources in its field of view, one filter at a time. The filters used in each epoch are reported in Table 2. Exposure times for MOSFIRE observations are given for each filter. The data reduction and spectral extraction were performed using MOSFIRE's data reduction pipeline,¹¹ and telluric and flux calibrations were performed using xtellcor. We caution that the accuracy of flux calibration between different epochs and different instrument can vary up to a factor of 2–3. Galaxy host contamination also varied among spectra due to the different slit widths and observing conditions. However, all subsequent analyses of the spectra will not rely on absolute flux calibration. All 1–2.5 μm spectra are shown in Figure 4.

Zoom In Zoom Out Reset image size

Finally, we obtained 3.0–3.9 μm spectroscopy using the Gemini Near-InfraRed Spectrograph (GNIRS) on the Gemini North telescope as part of a fast turnaround program (GN-2018A-FT-108, PI Tinyanont). Ground-based observations in this wavelength is difficult owing to the very strong thermal background from both the sky and the telescope. For these observations, we used the long red camera with the pixel scale of 005/pixel and the 015 slit with the 10 l/mm grating for the resolving power of R ∼ 1200. The wavelength range captured in our spectra was 3.1–3.9 μm. To handle high background, we set the detector bias to the deep-well mode to increase the saturation threshold and the read-out mode to match the read-out noise to the sky background noise. The observations were taken on 2018 August 16–17, with the total integration time of 23 and 68 minutes on the respective nights. We found, however, that only the first 17 images from the first night contain the SN trace. Each exposure was 30 s, 2 coadds, in ABBA dithering pattern along the slit. The data were reduced using Gemini's IRAF-based data reduction pipeline.¹² The telluric and flux calibrations were performed using xtellcor with HIP 114714 (HD 219290) as a standard star. We caution that spectral features present are likely spurious, and the spectra simply probed the continuum. We overplotted the binned spectrum in Figure 6, but did not use it to fit the SED.

Thermal emission arises from dust grains at different temperatures and different locations in the SN. To roughly determine this, we first fit a simple blackbody to the IR SED of SN 2014C at 1623 days, the epoch at which we have the 10 μm data and for which we will perform full SED fitting. We find that the near-IR (1–2.5 μm) flux is well explained by a hot component at T ∼ 1000 K located at r ∼ 10¹⁵ cm. The 3–5 μm flux arises from a warm component with T ∼ 500 K and r ∼ 4 × 10¹⁶ cm. Finally, the 10 μm flux can be explained by a cold component with T ∼ 300 K and r ∼ 1.5 × 10¹⁷ cm. We will argue later that the 10 μm flux is more likely from silicate dust emission feature and not a cold dust component. We note that the blackbody radii are the physical radii of optically thick dust clouds. In reality, the CSM dust is optically thin, and the blackbody radii are lower limits to the actual location of the dust grains. The shock front at this epoch, which is determined from the model that we will discuss in Section 3.4, is at r_sh = 1.2 × 10¹⁷ cm, suggesting that the hot component arises from likely newly formed dust inside the ejecta while the warm component is likely associated with pre-existing CSM dust heated by the shock either radiatively or collisionally.

We determine the composition of the dust in the CSM around SN 2014C by fitting its IR SED with different dust models. Following T16 and Fox et al. (2010), we fit the SED with a modified blackbody spectrum,

$$\begin{eqnarray}&&{F}_{\nu }=\displaystyle \frac{{M}_{\mathrm{dust}}{B}_{\nu }({T}_{\mathrm{dust}}){\kappa }_{\nu }(a)}{{d}^{2}},\end{eqnarray} \tag{ 1 }$$

where M_dust is the total dust mass, B_ν is the Planck function, T_dust is the dust temperature, and κ_ν(a) is the dust grain opacity, which is a function of dust grain size, a. This dust SED assumes that the emitting dust is optically thin and at a single temperature. We will check the optical depth assumption after calculating the dust mass. We note that this equation poses no geometric requirements on the distribution of the dust grains. The dust emissivity (Q) and consequently the opacity (κ) for carbonaceous and silicate grains are from Draine & Lee (1984) and Laor & Draine (1993), where we only consider graphite for carbonaceous grains. We convert the tabulated emissivity into opacity using the relation κ = 3Q/4ρ_bulka, where the bulk densities are ρ_bulk = 2.2 and 3 g cm⁻³ for the graphite and silicate grains, respectively. Amorphous carbon has higher emissivity, requiring less mass to explain the same IR flux, but the general shape of the SED is similar to that of graphite. We assume the grain size of 0.1 μm, noting that the emissivities from different grain sizes are degenerate in the IR, unless we include the unlikely large 1 μm grains (see Figure 4, top panel from Fox et al. 2010).

Figure 5 shows the SED of SN 2014C with near-IR photometry from P200/WIRC, mid-IR photometry from Spitzer/IRAC and Gemini/NIRI, and finally mid-IR (9.7 μm) photometry from Subaru/COMICS. All data were taken between 1622 and 1626 days post-maximum. The Spitzer data were interpolated from days 1515 and 1889 using a power law. The two panels of Figure 5 compare two different models, which provide different explanations for the 9.7 μm flux: copious cold carbonaceous dust (left) and the silicate dust's feature (right).

Zoom In Zoom Out Reset image size

We first consider a scenario where all CSM dust around SN 2014C is carbonaceous (Figure 5, left). This composition choice is based on prior observations, showing that no other interacting SNe with mid-IR observations showed silicate features. In this model, the 9.7 μm flux is dominated by cold carbonaceous dust at 276 K. We sample the parameter space with a Markov Chain Monte Carlo (MCMC) routine emcee (Foreman-Mackey et al. 2013) to derive the model parameters along with their robust uncertainties. More details on our MCMC sampling method are given in Appendix A.1, and the corner plot of this MCMC fitting is described in Figure 9. While this model fits our data well, the high mass of cold carbonaceous dust necessary to produce the observed 9.7 μm flux, 0.22 M_⊙, is problematic. First, 0.22 M_⊙ of carbonaceous dust at a minimum distance of r_bb = 1.5 × 10¹⁷ cm (the blackbody radius, derived earlier in this section for a cold dust component, also consistent with the shock radius derived later in Section 3.4) would have maximum optical depths, ${\tau }_{\nu }={M}_{d}{\kappa }_{\nu }/4\pi {r}_{\mathrm{bb}}^{2}=3$, 2, and 0.6 at 3.6, 4.5, and 9.7 μm respectively. This means that the dust could be not optically thin, resulting in even higher total dust mass than observed. Second, this dust mass is an order of magnitude higher than that dust mass observed in any SNe at this epoch (see, e.g., Figure 4 Gall et al. 2014) and is an order of magnitude higher than the dust mass inferred in SN 2006jd from the 3 to 10 μm observations at a similar epoch (Stritzinger et al. 2012). We also note that SN 2006jd had λ > 5 μm observations with neither signs of silicates nor cold (∼300 K) carbonaceous dust components that appear in SN 2014C. As mentioned earlier, the dust mass inferred here is the mass of the IR-emitting dust, which is only a lower limit to the total dust mass in the CSM. Assuming the canonical gas-to-dust mass ratio of 100, the lower limit of the emitting dust mass indicates that the mass of the gas shell containing the IR-emitting dust is ≳20 M_⊙, which is large in comparison to even the most massive Galactic LBV nebulae (Smith & Owocki 2006; Kochanek 2011; these authors assumed the same gas-to-dust ratio). Furthermore, we will show in Section 3.4 that the total CSM mass around SN 2014C could only be about 4–9 M_⊙. These are the major caveats to the purely carbonaceous dust scenario.

We now consider a second scenario, in which the 9.7 μm flux of SN 2014C is explained by the broad silicate feature around 10 μm (Figure 5, right). In this model, we fit the SED with two dust components. The first component is a single-temperature warm dust (T_warm) composed of both carbonaceous and silicate grains with a total mass of M_warm, a fraction f_Si of which is silicate. The second component is hot dust at the temperature T_hot composed of only carbonaceous dust with the total mass of M_hot. We used MCMC as described previously and in Appendix A.1 to derive the model parameters along with their robust uncertainties. The composition of this hot component does not matter because the flux at 10 μm from this component is negligible. The resulting parameters from the MCMC fits are ${T}_{\mathrm{warm}}={484}_{-13}^{+11}\,{\rm{K}}$, ${M}_{\mathrm{warm}}\,={5.2}_{-0.6}^{+0.7}\times {10}^{-3}\,{M}_{\odot }$, ${T}_{\mathrm{hot}}={1165}_{-100}^{+140}\,{\rm{K}}$, ${M}_{\mathrm{hot}}={4}_{-2}^{+4}\times {10}^{-6}\,{M}_{\odot }$, and ${f}_{\mathrm{Si}}={0.38}_{-0.08}^{+0.06}$. The corner plot of the MCMC fitting is presented in Figure 10. The uncertainties provided here only incorporate uncertainties in the measured flux and not in the dust emissivity. The near-IR photometry is well fitted by a hot, 1165 K, carbonaceous dust component. The total warm dust mass in this scenario, 5 × 10⁻³ M_⊙, is similar to that observed in SN 2006jd (Stritzinger et al. 2012). The dust cloud of this mass could also be optically thin if dust is located around the shock radius of 1.2 × 10¹⁷ cm; the maximum optical depths are 0.01, 0.01, and 0.1 at 3.6, 4.5, and 9.7 μm, respectively. Because of the more reasonable dust mass required to fit the SED, we favor this scenario in which the 9.7 μm is dominated by the broad silicate dust emission. Follow-up observations with more photometric bands around the 9.7 μm band are required to robustly identify the silicate feature. The mixture of silicate and carbonaceous dust grains required to explain the SED suggests that the CSM is inhomogeneous in its composition because a homogeneous medium with ${\rm{C}}/{\rm{O}}\ne 1$ will form only either carbonaceous or silicate dust. We will further discuss the implication on the origin of the CSM around SN 2014C in Section 4.

Figure 6 compares the IR SED of SN 2014C with the SEDs of all other interacting SNe for which observations beyond 5 μm are available in the literature. Apart from SN 2014C, only four other strongly interacting SNe IIn have been observed in the mid-IR. Three of these were observed at epochs comparable to the epoch at which SN 2014C was observed in the mid-IR. SN 2005ip was observed at 936 days post-explosion during the cryogenic phase of Spitzer. It was the only interacting SN with an InfraRed Spectrograph (IRS; Houck et al. 2004) mid-IR spectrum from 5 to 12 μm and IRAC photometry at 5.8 and 8 μm (Fox et al. 2010). SN 2006jd was observed at 1638 days, an epoch very similar to that of SN 2014C, with Spitzer/IRAC and WISE (Stritzinger et al. 2012). SN 2010jl was observed at 1279 days with Spitzer/IRAC and SOFIA/FORCAST (Herter et al. 2018) at 11.1 μm, resulting in a deep upper limit after 6400 s of total integration time (Williams & Fox 2015). In all three cases, the photometry and spectra from 1 to 10 μm are well fitted by purely carbonaceous dust models, which we overplot in Figure 6 using dust parameters from the literature. Lastly, SN 1995N was observed at more than 10 yr post-explosion by Spitzer/IRAC and WISE (Van Dyk 2013). Its SED shape differs markedly from those of the other three SNe observed at earlier epochs. The carbonaceous dust model that Van Dyk (2013) fitted to the data is shown in Figure 6. However, we note that the shallow slope from 3 to 10 μm cannot be fitted with a single-temperature dust model, regardless of composition, and would require a range of dust temperatures. Given the late epoch of the observation, one might also consider a nonthermal origin for the IR emission from SN 1995N, as the SED can also be described with a broken power law (F_ν ∝ ν⁻³; also overplotted) with a knee at around 12 μm. In addition to these H-rich interacting SNe, the H-poor interacting SN 2006jc (Ibn) was also observed beyond 5 μm with AKARI (Sakon et al. 2009). Its SED was also best fitted with a two-temperature amorphous carbon dust model. This comparison highlights SN 2014C's unique SED shape among other interacting SNe for which data are available beyond 5 μm at comparable epochs, showing for the first time an evidence for a silicate dust feature in the IR SED of an interacting SN. It also accentuates the need for observations of interacting SNe at late times, out to decades post-explosion, in the near- to mid-IR, which will be enabled by the upcoming James Webb Space Telescope (JWST).

Zoom In Zoom Out Reset image size

To track the evolution of the dust parameters, we fit the SED at other epochs with Spitzer/IRAC data, assuming a fixed ratio between carbonaceous and silicate dust. We do not fit for the hot, near-IR component, due to the lack of temporal coverage and its small contribution (<10%) to the IR luminosity. We run the MCMC fitting routine on all epochs and find that the posterior distributions of dust luminosity, temperature, and mass are Gaussian. Those parameters are shown in the top, middle, and bottom panels of Figure 7 and are listed in Table 3. At 1620 days, the epoch for which we have the full SED, we compare the values obtained for the warm dust component by fitting the full SED to those obtained from fitting only 3–5 μm data. This is to demonstrate that the discrepancy is minimal. In the top panel, we also show the full IR luminosity including both the hot and warm dust components to demonstrate that the hot component contributes minimally to the total IR luminosity. As a result, we use the luminosity of the warm dust component as an estimate of the bolometric luminosity of the SN.

Zoom In Zoom Out Reset image size

In this section, we fit the light curve of SN 2014C obtained in the last section with an analytic model from Moriya et al. (2013). The model assumes that the density profile in the SN ejecta follows a broken power law with ${\rho }_{\mathrm{ej}}\propto {r}^{-\delta ,-n}$ for the inner and outer ejecta, respectively. These ejecta interact with the CSM, whose density profile is ρ_CSM = Dr^−s. The interaction starts in the outer ejecta, and then transitions into the inner ejecta at time t_t post-explosion. Moriya et al. derived analytic expressions for the bolometric luminosity in the general CSM case where s and D are free parameters, for both t < t_t (Equation (23)) and t > t_t (Equation (26)). The latter scenario is applicable to our data taken after 600 days, when the flux in the Spitzer bands appear to dominate the SED. We fit the following equation (Equation (26) in Moriya et al. 2013) to our data:

$$\begin{eqnarray}&&L(t)=2\pi \epsilon {{Dr}}_{\mathrm{sh}}{\left(t\right)}^{2-s}{\left(\displaystyle \frac{(3-s){M}_{\mathrm{ej}}{\left(\tfrac{2{E}_{\mathrm{ej}}}{{M}_{\mathrm{ej}}}\right)}^{1/2}}{4\pi {{Dr}}_{\mathrm{sh}}{\left(t\right)}^{3-s}+(3-s){M}_{\mathrm{ej}}}\right)}^{3},\end{eqnarray} \tag{ 2 }$$

where the time dependence is in r_sh(t), which is numerically solved using their Equation (12):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{4\pi D}{4-s}{r}_{\mathrm{sh}}{\left(t\right)}^{4-s}+(3-s){M}_{\mathrm{ej}}{r}_{\mathrm{sh}}(t)\\ \quad -(3-s){M}_{\mathrm{ej}}{\left(\displaystyle \frac{2{E}_{\mathrm{ej}}}{{M}_{\mathrm{ej}}}\right)}^{1/2}t=0.\end{array}\end{eqnarray} \tag{ 3 }$$

We will later compute t_t to verify that our data satisfy t > t_t.

First, we perform an MCMC fit varying the CSM parameters: D and s from ρ_CSM = Dr^−s, and the explosion parameters M_ej and E_ej. We impose minimal physical constraints: M_ej, E_ej, D > 0, and the explosion velocity v_ej = (2E_ej/M_ej)^0.5 < c. The analytic model is only valid for 0 ≤ s < 3. We also try removing the degeneracy between M_ej and E_ej by imposing the observed photospheric velocity v_phot = (10E_ej/3M_ej)^0.5 = 13,000 km s⁻¹ (Margutti et al. 2017; Milisavljevic et al. 2015). In both cases, MCMC fits do not converge as the explosion parameters mainly affect the shape of the light curve at early times. The fitted explosion energy can grow arbitrarily large with a steep CSM density profile that makes very dense CSM close to the star.

Instead, we use the explosion parameters inferred from early-time observations: M_ej ∼ 1.7 M_⊙ and E_ej ∼ 1.8 × 10⁵¹ erg (Margutti et al. 2017), and only fit for the CSM parameters. The results we obtain, with 1σ uncertainty, are D = 10^14.9±0.2 (cgs) and s = 2.01 ± 0.01, which indicate that the CSM profile is wind driven. Recall that the parameter D sets the absolute density scale of the CSM. In the wind-driven scenario (s = 2), $D=\dot{M}/(4\pi {v}_{w})$, where $\dot{M}$ and v_w are the mass-loss rate and the wind velocity, respectively, and $\dot{M}/{v}_{w}$ is dependent on the explosion parameters. To derive the exact dependence, we solved Equation (35) and (36) in Moriya et al. (2013) to eliminate and obtained $D\,\propto \dot{M}/{v}_{w}\propto {M}_{\mathrm{ej}}^{3/2}\,{E}_{\mathrm{ej}}^{-1/2}$. From our fitted distribution for D and this relation, we can write the mass-loss rate as a function of v_w and the explosion parameters as

$$\begin{eqnarray}\begin{array}{rcl}\dot{M} & = & ({1.7}_{-0.6}^{+0.9})\times {10}^{-3}\,{M}_{\odot }{\mathrm{yr}}^{-1}\left(\displaystyle \frac{{v}_{w}}{100\,\mathrm{km}\,\ {{\rm{s}}}^{-1}}\right)\\ & & \times \,{\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{1.7{M}_{\odot }}\right)}^{1.5}{\left(\displaystyle \frac{{E}_{\mathrm{ej}}}{1.8\times {10}^{51}\mathrm{erg}}\right)}^{-0.5}.\end{array}\end{eqnarray} \tag{ 4 }$$

To check that we are in the t > t_t limit, we compute the transition time using Equation (34) in Moriya et al. (2013) and found that t_t = 345 ± 12 days. Hence, our observations from t > 600 days are well beyond this limit.

The density profile derived above allows us to estimate the total mass contained in the CSM. With the wind-driven density profile ρ ∝ r⁻², the gas mass contained up to radius r is linearly proportional to r. The total (gas) mass of the wind-driven component of the CSM contained up to radius r is simply M_CSM = 4πDr. Because D = 10^14.9±0.2 g cm⁻¹, the wind-driven portion of the CSM contains ${M}_{\mathrm{CSM}}={5}_{-2}^{+3}\times {10}^{-18}{\rm{\Delta }}r\,{M}_{\odot }\,{\mathrm{cm}}^{-1}$. As such, even if the outer CSM extends out to 10 times the current shock location, 10¹⁸ cm, the total wind-driven CSM mass is only about ${5}_{-2}^{+3}\,{M}_{\odot }$. More likely, the region of the CSM containing the dust we observe only extends out the order of its current size (∼2 × 10¹⁷ cm), with the total mass of about 1 M_⊙. Along with ≈1 M_⊙ in the dense CSM shell (Margutti et al. 2017), the total CSM mass around SN 2014C is about ${2}_{-0.4}^{+0.6}\,{M}_{\odot }$. This total CSM mass constraint presents a further issue for the pure carbonaceous dust scenario, which needs 0.22 M_⊙ of carbonaceous dust to explain the 9.7 μm flux observed at 1620 days. Recall that this dust mass is a lower limit, assuming that dust is optically thin and not clumpy, of the total mass of dust emitting in the IR at that epoch. That quantity is itself a lower limit of the total mass of dust in the CSM because not all dust is illuminated by the shock interaction at that epoch. The very conservative upper limit of gas-to-dust ratio is then ${9}_{-2}^{+3}:1$, which is at odds with dusty media observed around Galactic massive stars (e.g., 38 ± 15:1 in LBV WRAY 15-751; Vamvatira-Nakou et al. 2013) and SN remnants (e.g., the Crab Nebula has gas-to-dust ratio of 26–39:1; Owen & Barlow 2015).

Figure 8 (left) shows the evolution of the line profile of the broad He I 1.083 μm line. The 2.058 μm line has a similar profile at all epochs, but the S/N is low as it is a much weaker line. We analyze spectra from days 1354 and 1693 because they have the highest S/N, and the other epochs that cover the Y band (1305 and 1621 days) are close in time to these two epochs. We model the profile of the 1.083 μm line with a combination of Gaussian components and find that both epochs are well fit by four Gaussian components. Figure 8 (middle and right) shows the fitting results. There is a broad component approximately at rest, marked as "a," with a FWHM velocity of ≈2000 km s⁻¹. This component is likely from the He-rich ejecta of the SN. Second, there are two intermediate components, marked "b" and "c," at the central velocities of −4000 and about 0 km s⁻¹. They are both with a common FWHM of 500 km s⁻¹, suggesting that these intermediate lines are from the shocked CSM. Lastly, one narrow component (n1) with FWHM < 100 km s⁻¹ is present. The n2 component is the hydrogen 1.094 μm line. The FWHM of all resolved lines (not n1 and n2) decreases by about 10% between the two epochs. For comparison, the CSM interaction model used in Section 3.4 also tracks the shock velocity. In the wind-driven case, we can differentiate Equation (19) in Moriya et al. (2013) with respect to t and get the shock velocity

$$\begin{eqnarray}&&{v}_{\mathrm{sh}}(t)=\sqrt{\displaystyle \frac{2{E}_{\mathrm{ej}}/{M}_{\mathrm{ej}}}{1+2{at}}}\quad \mathrm{where}\quad a=\sqrt{\displaystyle \frac{2{E}_{\mathrm{ej}}}{{M}_{\mathrm{ej}}^{3}}}\displaystyle \frac{\dot{M}}{{v}_{w}}.\end{eqnarray} \tag{ 5 }$$

According to this equation, the shock velocity would decrease by 5% between 1354 and 1693 days, given the assumed ejecta mass, energy, and the fitted mass-loss rate and wind velocity inferred from the light curve. This is roughly consistent with the deceleration observed in the ejecta. The "b" component of the spectral profile may arise from the same hot spot found in the VLBI imaging (Bietenholz et al. 2018).

Zoom In Zoom Out Reset image size

The CSM observed in SN 2014C clearly has several distinct components, both in density and composition. As shown already in the literature (Milisavljevic et al. 2015; Anderson et al. 2017; Margutti et al. 2017), the CSM has at least two density components: a low-density bubble up to ∼10¹⁶ cm, and a dense, hydrogen-rich shell starting at ∼5 × 10¹⁶ cm. Continued IR observations, presented in this work, show that the CSM extends out to at least 1.4 × 10¹⁷ cm based on the current shock location at the latest Spitzer epoch at 1920 days. The light-curve shape, fitted with a CSM interaction model, further shows that the density profile of the outer part of the CSM differs from the constant-density shell observed at early times. The outer CSM density profile is well explained by a steady wind-driven scenario (ρ_CSM ∝ r⁻²) with a high mass-loss rate of $\dot{M}\approx {10}^{-3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The evolution of the observed dust mass also shows two distinct density profiles in the CSM. Before 750 days post-explosion, the observed dust mass is growing as M_dust ∝ t², indicating a CSM shell with constant density. After 750 days, the observed dust mass remains constant, indicating the CSM with a wind-driven (ρ_CSM ∝ r⁻²) density profile. In addition to the distinct density components, we infer from the SED fitting that the CSM likely contains both carbonaceous and silicate grains, which form in different environments (C/O ratio > 1 and <1, respectively), suggesting that the CSM was formed by stellar outflows with different compositions.

One possible explanation for the density and chemical inhomogeneity is that the CSM is created by a progenitor that is transitioning from the oxygen-rich phase (e.g., RSG or LBV) to the carbon-rich (e.g., WC) phase, similar to what has been suggested by Milisavljevic et al. (2015). In this scenario, the dense shell is formed by the fast WR wind sweeping up the slower RSG (or LBV) wind, which can mix in carbonaceous dust. There are three caveats of this scenario. (i) The WR phase would only last for ∼10 yr because the wind velocity is ∼1000 km s⁻¹ and the CSM shell is at ∼10¹⁶ cm. This is much shorter than the expected lifetime of a WR. Furthermore, LBVs (or RSG phase) are not expected to transition directly into the WC phase, but have to go through the WN phase first, making this scenario even less likely. This caveat was mentioned in Milisavljevic et al. (2015). (ii) A single WR progenitor is incompatible with the low ejecta mass of ∼1.7 M_⊙ inferred by Margutti et al. (2017). (iii) The wind parameter of $\dot{M}/{v}_{w}=1.7\pm 0.1\times {10}^{-3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}/100\,\mathrm{km}\,{{\rm{s}}}^{-1}$, inferred from our light-curve fitting, is inconsistent with the LBV quiescent wind phase. It is more consistent with either the high end of extreme RSG mass loss or the binary Roche-lobe overflow (RLOF) scenario (Smith 2017).

These caveats can be resolved by considering a binary progenitor system to SN 2014C. In this scenario, the H-rich envelope of the progenitor star is stripped by binary interaction, not by its line-driven stellar wind (e.g., Podsiadlowski et al. 1992). As a result, the progenitor star need not be massive to be able to shed its hydrogen envelope, which is consistent with the ejected mass inferred from observations. Margutti et al. (2017) discussed this scenario in detail as one of the possible origins of the dense CSM shell around SN 2014C (another scenario being late-stage nuclear-burning instability). They found from binary population modeling that 3%–10% of SNe Ib/c progenitor went through a common-envelope mass-loss scenario in which the stripped H-rich envelope can remain around the progenitor system long enough to interact with the SN. Such a binary progenitor system can also explain the chemical inhomogeneity and the density profile in the outer CSM. The density profile created by a binary system, assuming no interactions between two winds, is ${\rho }_{\mathrm{CSM}}(r)=({\dot{M}}_{1}/{v}_{1}+{\dot{M}}_{2}/{v}_{2})/(4\pi {r}^{2})$, where ${\dot{M}}_{i}$ and v_i are the mass-loss rate and wind velocity of the two components, and r is much larger than the binary separation. This is consistent with the observed density profile presented in this work. This scenario is also in line with the growing consensus that most stripped-envelope SNe are created via binary interaction, as opposed to single star mass loss (e.g., Smith 2014; De Marco & Izzard 2017). SN 2014C stands out from the rest of SNe Ib/c, perhaps due to the progenitor's binary parameter that caused the mass-ejection event to occur right before core collapse of one of the components.

There is a potential Galactic analog to the progenitor system of SN 2014C: the massive binary system RY Scuti (Hjellming et al. 1973; Gehrz et al. 1995; Smith et al. 1999; Gehrz et al. 2001; Smith et al. 2002, 2011a). RY Scuti is an interacting binary caught shortly after major mass-transfer events. It consists of a primary, mass donor, O9/B0 supergiant star with a current mass of 8 M_⊙ and a secondary, mass gainer, O5 star with the current mass of 30 M_⊙. The initial masses for the primary and secondary stars were estimated to be about 20–25 and 15–20 M_⊙. The He core mass of the primary star would have been roughly 6–8 M_⊙ (Woosley & Weaver 1986), indicating that it had lost most of its hydrogen envelope and would soon become a WR star. High-resolution imaging in the optical with Hubble Space Telescope (HST) and in the infrared using adaptive optics (AO) revealed detached, expanding tori of CSM ejected in separate eruptive mass-loss events (Smith et al. 1999, 2002, 2011a). The inner gas-rich torus was 3.4 × 10¹⁶ cm from the star and was detected primarily in the optical. The outer dusty torus was 5.4 × 10¹⁶ cm from the star and was detected primarily in the IR. Proper motion measurements from multiepoch observations 6 (12) yr apart in the infrared (optical) showed that the radial expansion velocity of the gas and dust tori was 40–50 km s⁻¹. The two tori had different velocities and were ejected separately, about 100 yr apart (Smith et al. 2011a).

The SN from the primary star of RY Scuti may resemble SN 2014C in the following aspects. First, the SN will initially be either of Type IIb or Ib/c, depending on how much hydrogen envelope is left. The SN will transition into a SN 2014C-like event once the SN shock reaches the location of the dusty tori. The strength of the interaction, however, may be much lower depending on the mass in the CSM. The outer dusty torus around RY Scuti was 5.4 × 10¹⁶ cm away from the star (measured from imaging), which is comparable to the location of the dense CSM shell around SN 2014C (Margutti et al. 2017). Hence, the timescale at which the shock interaction commences would be similar to that observed in SN 2014C (about 100 days post-explosion). Second, the toroidal geometry of the CSM will allow the SN shock to propagate unimpeded, explaining the constant shock velocity despite CSM interaction (Bietenholz et al. 2018). Third, the IR SED of the outer dusty torus observed by the Keck Long Wave Spectrometer (LWS) showed a strong silicate dust feature at 10 μm, similar to what we inferred in SN 2014C (Gehrz et al. 2001). Moreover, Gehrz et al. (2001) invoked a dust model with a mix of carbonaceous and astronomical silicate dust grains to explain the SED shape of RY Scuti, similar to what we did for SN 2014C. The temperature of the dust grains was slightly lower to what was observed in SN 2014C (300–400 K as opposed to 500 K), which is not unexpected given the lower luminosity from RY Scuti (≈2 × 10⁶ L_⊙) as opposed to what would be expected in a SN explosion. While RY Scuti and SN 2014C's progenitor may be similar in many aspects, they differ sharply in the CSM mass. SN 2014C has of order 1 M_⊙ in the CSM with about 5 × 10⁻³ M_⊙ of dust; meanwhile, RY Scuti has only 1.4 × 10⁻⁶ M_⊙ of dust (Gehrz et al. 2001). This led Smith et al. (2011a) to conclude that the SN from RY Scuti's primary star will be an ordinary SN Ib/c and not an interacting SN IIn. They noted, however, that RY Scuti's mass transfer was very conservative, because about 15 M_⊙ has been transferred to the secondary star with only about 0.003 M_⊙ lost into the CSM (gas mass of the inner torus; Smith et al. 2002). The progenitor system to SN 2014C may be a binary system much like RY Scuti, but with mass ratio and separation such that the mass transfer was much less conservative, resulting in much more mass being lost into the CSM.

SN 2014C's prolonged IR emission reveals two important characteristics of the CSM. First, it likely contains a mixture of carbonaceous and silicate dust inferred from the SED including 9.7 μm photometry. If confirmed, this is the first identification of silicate dust in an interacting SN. Alternatively, the SED could also be fitted with purely carbonaceous dust with a cold, 300 K, component containing 0.22 M_⊙ of dust grains. That scenario would be unprecedented as well. Second, the light-curve fitting shows that the CSM extends out to at least 1.4 × 10¹⁷ cm, with the wind-driven density profile in the outer part. This component is in addition to the inner low-density bubble and the dense shell with a constant density identified earlier in the literature. The chemical inhomogeneity of the CSM suggests that the progenitor system of SN 2014C is likely a binary with one component producing an oxygen-rich outflow and another component carbon-rich outflow. While most SNe Ib/c are not observed to interact with their lost hydrogen envelope, SN 2014C may present the extreme end of the population, where the binary-driven mass-ejection event that strips the progenitor's envelope happens just before the core collapse, leaving a dense, detached H-rich CSM to interact with the SN shock. The progenitor system that produced SN 2014C may look very similar to RY Scuti, but with much more mass ejected from the binary system into the CSM. This work underscores the importance of late-time IR observations to constrain CSM properties of interacting SNe, as the outer part of the CSM can only be probed by late-time interactions, which emit mostly in the IR. Assuming that this CSM does not abruptly end at this radius, the interaction will still be ongoing by the time that JWST becomes operational. It will enable IR spectroscopic observations of this SN and of other interacting SNe, which will help us understand SN 2014C's place in the continuum of SN types from Ib/c to IIn, and to construct a complete picture of how mass loss operates at the end of the life of a massive star.

To determine best-fit values for the dust parameters in the SED fittings in Section 3.1 and the CSM parameters in the light-curve fitting in Section 3.4, we used the emcee package to run MCMC on our data. For the first model with all carbonaceous dust, there were six parameters: M_{cold,warm,hot} and T_{cold,warm,hot}. The dust composition was purely carbonaceous, and the grain size was 0.1 μm. The initial values of these parameters were obtained by running a more traditional least-squares fit on the data using the curve_fit routine in the scipy.optimize package. The prior distribution for all parameters is uniform, with only positive values allowed for mass and temperatures. We ran MCMC using 400 walkers for 1000 steps and determined that the convergence happened after 200 steps. Figure 9 shows the resulting posterior distribution of the six parameters, along with the median and ±1σ values.

Zoom In Zoom Out Reset image size

For the second model with a mixture of carbonaceous and silicate dust, there were five free parameters: M_warm, T_warm, M_hot, T_hot, and f_Si. The prior distribution was also uniform with positive mass and temperatures. The silicate fraction f_Si could be between 0 and 1. The initial values were also obtained with a least-squares fit using curve_fit. We ran MCMC using 400 walkers for 2000 steps, and determined that the convergence happened after 1000 steps. We note that the temperature and mass of the hot component are not as well constrained because the near-IR SED deviates from the dust model, potentially due to poor background subtraction. Figure 10 shows the resulting posterior distribution of the five parameters, along with the median and ±1σ values.

Zoom In Zoom Out Reset image size

As discussed in the main text, we found that the MCMC fit can only constrain CSM density parameters, namely D and s from ρ_CSM = Dr^−s, and not the explosion parameters (E_ej and M_ej). We considered a model with D and s as free parameters. Similarly to the SED fit, the initial values were determined using least-squares fitting. We also used a uniform prior distribution with physical parameters set to positive numbers, and 0 ≤ s < 3. Figure 11 shows the corner plots for the fit.

Zoom In Zoom Out Reset image size