Supernova 2018cuf: A Type IIP Supernova with a Slow Fall from Plateau

The Milky Way line-of-sight reddening toward SN 2018cuf is E(B–V)_MW = 0.0273 ± 0.0003 mag (Schlafly & Finkbeiner 2011). This low extinction value is also supported by the lack of NaID lines from the Milky Way in the Magellan/MIKE spectrum taken on 2018 July 12 and shown in Figure 3. The equivalent width (EW) of the NaID line is often used to estimate the SN reddening with the assumption that it is a good tracer of gas and dust (Munari & Zwitter 1997; Poznanski et al. 2012). The measured EWs of the host galaxy NaID λ5890 (D₂) and NaID λ5896 (D₁) are 0.677 Å and 0.649 Å, respectively. The intensity ratio of D₂ to D₁ (D₂/D₁ ∼ 1) is far from the typical value of 2 we usually observe (Munari & Zwitter 1997), suggesting that at least D₂ may be saturated (see Figure 3). Using only D₁, we found a host galaxy extinction of E(B–V)_host = 0.699 ± 0.17 mag.

Phillips et al. (2013) suggested that the most accurate predictor of extinction is the diffuse interstellar band absorption feature at 5780 Å. However, this feature is not clearly present in our high-resolution spectrum of SN 2018cuf, suggesting the host galaxy extinction is low, which is inconsistent with the high host reddening derived from NaID lines. Munari & Zwitter (1997) found that [K i] λ7699 can be a better reddening indicator if NaID lines are saturated, so we decided to use this line to estimate the reddening from the host galaxy. The EW of [K i] λ7699 was measured to be 0.03 Å, which corresponds to a host galaxy extinction of E(B–V)_host = 0.11 ± 0.01 mag (Munari & Zwitter 1997). As a sanity check, we also compared the dereddened B − V color evolution of SN 2018cuf to a sample of similar SNe II with published reddening estimates. This includes SN 1993A (Anderson et al. 2014; Galbany et al. 2016b), SN 1999gi (Leonard et al. 2002b), SN 2003iq (Faran et al. 2014), SN 2003bn (Anderson et al. 2014; Galbany et al. 2016b), SN 2003ef (Anderson et al. 2014; Galbany et al. 2016b), SN 2003T (Anderson et al. 2014; Galbany et al. 2016b), SN 2009ib (Takáts et al. 2015), and SN 2012A (Tomasella et al. 2013), as is shown in Figure 6. SN 2018cuf has a similar V-band light-curve slope after maximum with these selected SNe. de Jaeger et al. (2018) found that the color evolution of SNe is related to the slope of the V-band light curve, so these selected SNe should have consistent colors with SN 2018cuf after dereddening. We found that an E(B–V)_host ≈0.11 mag gives us a consistent color evolution with the other objects. Therefore throughout this paper we will adopt an E(B–V)_tot = 0.1373 ± 0.0103 mag, as well as an R_V = 3.1 (Cardelli et al. 1989).

Zoom In Zoom Out Reset image size

The disagreement between host reddening values obtained from the NaID lines and direct color comparisons to similar objects is not a unique problem. Leonard et al. (2002a) found a similar situation for SN 1999em, i.e., the EW of the sodium lines suggested a high reddening for SN 1999em, but a low value was assumed based on color comparisons. Phillips et al. (2013) also found that NaID gives unreasonably high reddening for some of their objects, while the KI line gives a reddening that is consistent with the reddening derived from SN colors.

MUSE (Bacon et al. 2010) integral field unit observations of IC 5092 were taken on 2019 April 12, as a part of the All-weather MUSE Supernova Integral-field Nearby Galaxies (Galbany et al. 2016a) survey. MUSE is mounted to the 8.2 m Yepun UT4 Very Large Telescope, with a field of view of 1' × 1' and 02 × 02 spatial elements, small enough to sample the PSF. See Figure 1 for an outline of the MUSE footprint. The spectral coverage is from 4750 to 9300 Å, with a spectral resolution that ranges from R ≃ 3500 in the blue end to ≃1700 in the red end. Four 580 s exposures (2320 s total exposure time), rotating 90° between frames, were taken centered on the southwest side of the galaxy, which covered the SN position and its environment.

We extracted a 36 aperture spectrum centered at the SN position (corresponding to an ∼800 pc diameter) to study the properties of the environment. The resulting spectrum is shown in Figure 4. MUSE observations were performed 293 days after SN 2018cuf's explosion, and some SN features were still visible in the spectrum, with the most pronounced being a broad Balmer Hα emission, in addition to an H II region spectrum with narrow emission lines. To measure the flux of the strongest ionized gas emission lines in that region ([N ii] λ6548, Hα, and [N ii] λ6583), we excluded the SN broad component by fitting four Gaussians, three narrow and one broad, simultaneously. The bluer region of the spectrum was not strongly contaminated by SN features, and we fit single Gaussians to measure the narrow Hβ and [O iii] λ5007 emission line fluxes from the ionized gas.

An estimate of the reddening can be obtained from the line-of-sight gas column by the ratio of the Balmer lines, assuming a case B recombination (Osterbrock & Ferland 2006) and a theoretical ratio of Hα/Hβ = 2.86. Our lines present a ratio of 4.54, which corresponds to E(B–V) = 0.399 ± 0.021 mag. This value is not consistent with the reddening estimated from our color comparison (Figure 6) and would make the light curves of SN 2018cuf significantly bluer and brighter than those of similar SNe II. A possible explanation for this disagreement is that the SN is in front of the H II region and not influenced by the dust, but the MUSE measurement gets the full column of gas.

With the host galaxy reddening-corrected fluxes we estimated the SN environmental oxygen abundance (O/H) by using the N2 and O3N2 calibrators from Pettini & Pagel (2004). We obtained a consistent oxygen abundance of 12+log(O/H) = 8.71 ± 0.07 dex and 12+log(O/H) = 8.72 ± 0.08 dex with the N2 and O3N2 calibrators, respectively, both consistent with solar abundance (Asplund et al. 2009). We used the Hα luminosity to estimate the star formation rate (SFR) at the SN location using the expression provided by Kennicutt (1998). We obtained an SFR of 0.0014 ± 0.0001 M_⊙ yr⁻¹ and an SFR intensity of 0.0027 ± 0.0001 ${M}_{\odot }\ {\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$. To understand where SN 2018cuf stands in the SNe II group, we compared the values we derived above with the host properties of all SNe II in the PMAS/PPak Integral-field Supernova Hosts Compilation (PISCO) sample (Galbany et al. 2018).³² The average host oxygen abundance and SFR intensity for all PISCO Type II hosts are 12+log(O/H) = 8.53 ± 0.062 dex and 0.013 ± 0.0014 ${M}_{\odot }\ {\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$, respectively, suggesting that the region around SN 2018cuf has a higher oxygen abundance but a lower SFR intensity than the average of SNe II.

The distance to IC 5092 is not well constrained since it has only been measured using the Tully–Fisher relation (Mathewson et al. 1992; Willick et al. 1997) to be 32.0 ± 5.8 Mpc. While the Tully–Fisher relation can be used to measure distances to most spiral galaxies, the intrinsic scatter hinders the accuracy of the measurement for a single galaxy (Czerny et al. 2018). One commonly used approach to independently measure distances to SNe II is the expanding photosphere method (EPM), although it requires the object to have well-sampled light curves and spectra. The EPM was first developed by Kirshner & Kwan (1974) to calculate the distance to SNe IIP based on the Baade (1926) method. Assuming that the photosphere is expanding freely and spherically, we can obtain the distance from the linear relation between the angular radius and the expanding velocity of the photosphere using the function

$$\begin{eqnarray}&&t=D\left(\displaystyle \frac{\theta }{{v}_{\mathrm{phot}}}\right)+{t}_{0}\end{eqnarray} \tag{ 1 }$$

where D is the distance, t₀ is the explosion epoch, θ is the radius of the photosphere (in angular units), and v_phot is the velocity of the photosphere. Assuming that the photosphere radiates as a dilute blackbody, we combined the multiband photometry to simultaneously derive the angular size (θ) and color temperature (T_c) by minimizing the equation

$$\begin{eqnarray}&&\epsilon =\displaystyle \sum _{\nu \in S}{\{{m}_{\nu }+5\mathrm{log}[\theta \xi ({T}_{{\rm{c}}})]-{A}_{\nu }-{b}_{\nu }({T}_{{\rm{c}}})\}}^{2}\end{eqnarray} \tag{ 2 }$$

where ξ and b_ν are the dilution factor and synthetic magnitude, respectively, and both of them can be treated as a function of T_c (Hamuy et al. 2001; Dessart & Hillier 2005), A_ν is the reddening, m_ν is the observed magnitude, and S is the filter subsets, i.e., {BV}, {BVI}, and {VI}. We estimated the photospheric velocity by measuring the minimum of the [Fe ii] λ5169 P Cygni profile. To accurately estimate the error on this measurement and avoid noise-induced local minima, we smoothed the spectra with Savitzky–Golay filters (Savitzky & Golay 1964; Poznanski et al. 2010) with different widths, deriving the photospheric velocity for each width. For our distance measurement we used the mean and standard deviation of these velocity measurements. After ∼40 days, the relation between $\theta /v$ and t was clearly nonlinear (Jones et al. 2009), and for this reason we only used four early spectra with clear [Fe ii] λ5169 detection and interpolated the photometry data to the corresponding epoch. The measured velocities are listed in Table 1. In order to use the dilution factor derived by Dessart & Hillier (2005), we converted the rp and ip magnitude to I magnitudes by using the equations given by Lupton et al. (2005). The results for the three filter subsets {BV}, {BVI}, and {VI} are presented in Figure 7. From these measurements, we obtained distances of 43.3 ± 5.8 Mpc, 40.6 ± 5.2 Mpc, and 41.3 ± 7.4 Mpc, respectively, and the weighted average was calculated to be 41.8 ± 5.8 Mpc by using the method described in Schmelling (1995). In the rest of the paper, we will adopt this value for the analysis.

Zoom In Zoom Out Reset image size

We derived the explosion epoch from the EPM analysis, obtaining similar values from each of the three filter subsets used: JD 2,458,293.04 ± 2.88 days in {BV}, JD 2,458,292.95 ± 2.82 days in {BVI}, and JD 2,458,292.02 ± 4.22 days in {VI}. The weighted average of these measurements is JD 2,458,292.81 ± 3.08, which we adopt as the explosion epoch throughout this paper. We note that this is consistent with the tight constraints of the DLT40 survey, which place the explosion epoch between JD 2,458,291.74456 (the last nondetection) and JD 2,458,292.8609 (the first detection, which is just 0.05 d after the estimated explosion epoch).

As an independent check, we also estimated the explosion epoch by matching the spectra of SN 2018cuf with the spectral templates in the Supernova Identification (SNID) code (Blondin & Tonry 2007). This method has been used by Anderson et al. (2014) and Gutiérrez et al. (2017) to constrain the explosion epochs of a sample of SNe II. Gutiérrez et al. (2017) found that with the addition of new spectral templates to the SNID database, the explosion epoch derived from spectral matching may constrain the explosion to within 3.9 days. Following the work of Gutiérrez et al. (2017), we fixed the fitting range in SNID to 3500–6000 Å since the blue end of the spectrum contains more information about the SN and evolves more consistently with time for SNe II. Fixing the explosion epoch to JD 2,458,291.91 (from the EPM), we compared the spectra at 12.10, 19.06, 22.99, and 27.29 days with the SNID templates, where the explosion epochs are given by the EPM. The top five matches were then averaged to compute the epoch of the spectra, and the error was given by the standard deviation. The epochs of the spectra derived from this method are 10.84 ± 1.87, 17.82 ± 4.74, 25.02 ± 4.65, and 31.74 ± 6.79 days, respectively, consistent with the spectral epochs inferred from the EPM.

The optical spectroscopic evolution of SN 2018cuf is shown in Figure 4. The early spectrum shows a blue continuum with a broad Hα line clearly detected. Over time, the spectra become redder and develop hydrogen Balmer lines with P Cygni features. The [Fe ii] $\lambda \lambda \lambda$ 4924, 5018, and 5169 lines, good tracers for the photospheric velocity, can be seen after day 17. Other typical features such as [Ca ii] λλ3934, 3968, the Ca ii infrared triplet $\lambda \lambda \lambda$ 8498, 8542, 8662, and NaID λλ5890, 5896 also appear in emission as the SN evolves. During the nebular phase, strong [Ca ii] λλ7291, 7324 emission lines emerge along with [Fe ii] λ7155, [He i] λ6678, and [O i] λλ 6300, 6364.

Interestingly, from day 105 to day 174, a small notch appears on the Hα profile with a velocity of ∼1000 km s⁻¹, and its origin is unclear. One possibility is that this feature is from dust formation either in the ejecta or in circumstellar material (CSM) interaction. The signatures of dust formation have been detected in many Type IIn SNe. Type IIn SN2010jl shows notches or double-peaked profiles at an earlier stage and later shows more dominant blue-wings (see Extended Data Figure 3 of Gall et al. 2014). Type IIn SN1998S also shows a notch feature in its broad emission lines (Mauerhan & Smith 2012). However, for SN 2018cuf, this notch feature emerges starting at day 106, and the temperature of the ejecta may still be too hot for dust formation. On the other hand, the feature is not detected in the spectra after day 293, which is hard to reconcile with dust formation. By comparing Hα with other hydrogen line profiles, we do not find any evidence of red-side attenuation for lines that occur at bluer wavelengths, as is expected for dust creation. For these reasons, we cannot unambiguously attribute this feature to dust formation and equally rule out the possibility of dust formation.

The evolution of Hα and Hβ lines during the photospheric phase is shown in Figure 11. Starting at day 22, an extra absorption line can be seen on the blue side of the Hα and Hβ P Cygni absorption lines, becoming more obvious by day 34. These lines have been studied in many SNe II and have been most often associated with [Si ii] λ6355 when seen at early phases (<30 days) and with HV hydrogen when seen at a later phase (50–100 days; e.g., Chugai et al. 2007; SN 2005cs, Pastorello et al. 2006; SN 2009bw, Inserra et al. 2012; SN 2013ej, Valenti et al. 2014). This ambiguous absorption feature is often referred to as the "Cachito" feature (Gutiérrez et al. 2017). In the case of SN 2018cuf, because this feature appears at roughly 30 days, it is likely associated with HV hydrogen, an interpretation that is confirmed by the additional presence of the HV feature in Hβ at similar velocities.

Zoom In Zoom Out Reset image size

The NIR spectra from day 4 to day 94 are plotted in Figure 5 and show an evolution typical of SN II. The first spectrum at day 4 is nearly featureless with weak Paschen lines, but by day 15 these features have strengthened and [He i] λ10830 has also appeared. Both Pa α and Br γ lines can be seen in our spectra after day 58.

In general, the line evolution in NIR spectra is consistent across SNe II. However, Davis et al. (2019) pointed out that SNe II can be classified as spectroscopically strong or weak based on the pseudo-equivalent width (pEW) of the [He i] λ10830 absorption line and the features seen in the spectra. They found that SNe with weak He i (pEW < 50 Å) are slow-declining SNe IIP and that SNe with strong He I (pEW > 50 Å) correspond to the faster-declining Type IIL SNe class. Interestingly, SN 2018cuf seems to be an exception to this rule. The pEWs of He i absorption for days 39, 58, 81, and 94 are 4.9 Å, 25.2 Å, 82.5 Å, and 105 Å, respectively, which makes it hard to classify it as either a strong or weak SN based on the pEW alone. In addition to smaller pEWs, weak SNe usually show the Pγ/Sr II absorption feature at earlier epochs (∼20 days after explosion) and are accompanied by an HV He i feature. For SN 2018cuf, the Pγ/Sr II absorption feature shows up at day 18, consistent with a weak SN. Additionally, there is clearly an extra absorption feature on the blue side of He i (see the right panel of Figure 5 or Figure 11). Other than HV He I, this feature could also be explained as [C i] $\lambda 10693$. However, the lack of other C i lines in the NIR spectra makes it unlikely that this feature originates from [C i] $\lambda 10693$. We also note that the velocity of HV He I matches the velocity of HV Hα and HV Hβ in optical spectra, which further strengthens our conclusion that this feature can be interpreted as HV He I. Although the pEW of SN 2018cuf is greater than the 50 Å limit used in Davis et al. (2019), the presence of Pγ/Sr II at an early phase and of HV He I suggests that our object still falls into the weak SN II category. This implies that the 50 Å limit from Davis et al. (2019) is probably too low.

Chugai et al. (2007) proposed that HV absorption features, like those seen in SN 2018cuf, come from the interaction between the circumstellar wind and the SN ejecta. They argued that there are two physical origins of HV absorption: enhanced excitation of the outer layers of unshocked ejecta, which contributes to the shallow HV absorption in the blue side of Hα and [He i] λ10830, and the cold dense shell (CDS), which is responsible for the HV notch in the blue wing of Hα and Hβ. For the former case, the Hβ HV is not expected to be seen due to the low optical depth in Hβ line-forming regions, whereas in the latter case an HV Cachito can form in both Hα and Hβ. For SN 2018cuf, the presence of Cachito features in Hα, Hβ, and [He i] λ10830 supports the CDS interpretation but does not completely rule out the first scenario.

Paβ and Paγ were also investigated to look for HV features. However, the existence of other strong lines around Paγ makes it difficult to identify an HV feature if present, and there is no HV feature in the blue side of Paβ. Chevalier & Fransson (1994) suggested that the temperature of the CDS should be low enough that this region is dominated by low-ionization lines, which causes Paβ absorption to form in a low optical depth region and may explain the absence of HV features in Paβ.

The nebular phase of SNe II is driven by radioactive decay ⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe. If the γ-rays produced by this process are completely trapped by the ejecta, the bolometric luminosity at late times can be used to estimate the amount of ⁵⁶Ni. Since our photometry after ∼100 days does not cover the full spectral energy distribution (SED), we used two different methods to derive the ⁵⁶Ni mass. The first method is to calculate the pseudo-bolometric luminosity of SN 2018cuf and compare it with the pseudo-bolometric light curve of SN 1987A. Assuming SN 2018cuf and SN 1987A have the same normalized SED, the nickel mass is given by Spiro et al. (2014):

$$\begin{eqnarray}&&{M}_{\mathrm{Ni}}=0.075\times \displaystyle \frac{{L}_{\mathrm{SN}}}{{L}_{87A}}{M}_{\odot }\end{eqnarray} \tag{ 4 }$$

where L_SN and L_87A are the pseudobolometric luminosity of SN 2018cuf and SN 1987A, respectively. For the pseudo-bolometric luminosity, we followed the method described by Valenti et al. (2008). The observed magnitudes were converted to flux at each band and integrated by using Simpson's rule, which uses a quadratic polynomial to approximate the integral. The photometric data from day 135 to day 170 were used to calculate the pseudo-bolometric luminosities of SN 2018cuf and SN 1987A by using passbands {BVgri} and {BVRI}, respectively, resulting in a ⁵⁶Ni mass of ${0.037}_{-0.002}^{+0.003}{M}_{\odot }$.

An alternative approach to estimate the nickel mass is to compute a full-band bolometric light curve by performing a blackbody fit to all available filters at each photometric epoch and integrating the blackbody. The advantage of this method is that it does not require the assumption that SN 2018cuf and SN 1987A have the same normalized SED, although the approximation to a blackbody may not be completely valid due to the line blanketing in the UV bands. The ⁵⁶Ni mass derived from this approach is ${0.042}_{-0.008}^{+0.045}{M}_{\odot }$. Given the limitations of each method, we chose to use the pseudobolometric luminosity method to estimate the ⁵⁶Ni mass but took the difference between the results from the two methods as an indicator of the uncertainty of the measurement. The final nickel mass was conservatively estimated to be ${M}_{\mathrm{Ni}}=0.04\ (0.01)\,{M}_{\odot }$. By comparing the pseudobolometric light curve of SN 2018cuf with that of SN 1987A, we found that the decline rate of SN 2018cuf in the radioactive tail is either consistent with or slightly faster than ⁵⁶Co decay. It is hard to be sure which one is the case here due to the lack of data in the radioactive tail. If the decline rate of SN 2018cuf is slightly faster than ⁵⁶Co decay, the ⁵⁶Ni mass we derived here could be treated as a lower limit.

Progenitor mass is a fundamental parameter of an SN, and it can be constrained by using multiple techniques. In this section, we derive the progenitor mass of SN 2018cuf from nebular spectra and hydrodynamic light-curve modeling.

7.2.1. From Nebular Spectroscopy

During the nebular phase, spectra can provide useful information about the inner structure of an SN. At this stage, the ejecta has become optically thin, revealing the core nucleosynthesis products. The strength of the [O i] $\lambda \lambda 6300,6364$ doublet in the nebular spectra has been found to be a good indicator of progenitor mass (Jerkstrand et al. 2014). By comparing the intensities of [O i] $\lambda \lambda 6300,6364$ with theoretical models during this phase, the progenitor mass can be well constrained. Jerkstrand et al. (2014) modeled the nebular spectra for 12, 15, 19, and 25 M_⊙ progenitors at different phases. They started with the SN ejecta evolved and exploded with KEPLER (Woosley & Heger 2007) and used the spectral synthesis code described in Jerkstrand et al. (2011) to generate the model spectra.

Although we have six nebular spectra for SN 2018cuf taken from day 174 to day 483, four of them are contaminated by the host galaxy. Therefore, we only compare the nebular spectra of SN 2018cuf at day 174 and day 336 with the models computed by Jerkstrand et al. (2014) in Figure 12. We scaled the nebular spectra taken at day 174 and day 336 to the r-band photometry, and the models were scaled to the observed spectrum so that they had the same integrated flux. We found that the strength of O i in our spectrum is between the 12 M_⊙ and the 15 M_⊙ models, which implies the progenitor mass of SN 2018cuf is likely in this range.

Zoom In Zoom Out Reset image size

Synthetic nebular spectra can also be used to give an independent estimate of the nickel mass (Jerkstrand et al. 2018; Bostroem et al. 2019). By using the scale factors we used to scale the model spectra, the nickel mass can be derived using the following relation from Bostroem et al. (2019):

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{\mathrm{obs}}}{{F}_{\mathrm{mod}}}=\displaystyle \frac{{d}_{\mathrm{mod}}^{2}}{{d}_{\mathrm{obs}}^{2}}\displaystyle \frac{{M}_{{}^{56}{\mathrm{Ni}}_{\mathrm{obs}}}}{{M}_{{}^{56}{\mathrm{Ni}}_{\mathrm{mod}}}}\exp \left(\displaystyle \frac{{t}_{\mathrm{mod}}-{t}_{\mathrm{obs}}}{111.4}\right)\end{eqnarray} \tag{ 5 }$$

where F_obs is the total observed flux and F_mod is the total flux from the model spectrum. d_obs is the distance of the SN and d_mod = 5.5 Mpc is the distance used to compute the model; ${M}_{{}^{56}\mathrm{Ni}}$ indicates the nickel mass for the SN (${M}_{{}^{56}{\mathrm{Ni}}_{\mathrm{obs}}}$) and the model (${M}_{{}^{56}{\mathrm{Ni}}_{\mathrm{mod}}}=0.062{M}_{\odot }$), and t_obs and t_mod are the phase of the spectra for the observation and model, respectively. We then derived nickel masses of ${0.040}_{-0.003}^{+0.004}$ M_⊙ and ${0.060}_{-0.016}^{+0.006}$ M_⊙ for day 174 and day 336, respectively. These values are consistent with what we got in the previous subsection, where we measured the nickel mass from the radioactive decay tail photometry when the SN just falls from the plateau (days 135–170).

7.2.2. From Hydrodynamic Modeling

An alternative way of constraining the mass of the progenitor is to compare the light curves to hydrodynamic models (e.g., Utrobin & Chugai 2015, 2017; Morozova et al. 2017, 2018; Paxton et al. 2018; Goldberg et al. 2019; Martinez & Bersten 2019). We used the Supernova Explosion Code (SNEC; Morozova et al. 2015), an open-source hydrodynamic code for core-collapse SNe, to constrain the progenitor parameters of SN 2018cuf. SNEC assumes diffusive radiation transport and local thermodynamic equilibrium, which are valid assumptions from shock breakout through the end of the plateau. However, as the SN becomes nebular, this assumption breaks down. For this reason, we compared our light curve only out to ${t}_{\mathrm{PT}}=112.24$ days with the SNEC models. Our inputs of evolved progenitor stars for SNEC are the nonrotating solar metallicity RSG models generated from the kepler code and described in Sukhbold et al. (2016). A steady-state wind with a density profile

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{\dot{M}}{4\pi {r}^{2}{\nu }_{\mathrm{wind}}}=\displaystyle \frac{K}{{r}^{2}}\end{eqnarray} \tag{ 6 }$$

is also added above these models to explore the effect of CSM on light curves, where $\dot{M}$ is the wind mass-loss rate and ${\nu }_{\mathrm{wind}}$ is the wind velocity. We will use the parameter K to describe the constant wind density, which extends up to radius R_ext. For each explosion, SNEC takes a variety of progenitor and explosion parameters as input and then generates a bolometric light curve and, assuming blackbody radiation, the optical light curves. We followed the approach of Morozova et al. (2017), exploring variations in progenitor mass (M), nickel mass, explosion energy (E), K, and R_ext and fixing the nickel mass to ${M}_{\mathrm{Ni}}=0.04\ {M}_{\odot }$, which we obtained from the tail photometry. We note that the degree of ⁵⁶Ni mixing can also be a free parameter in SNEC models. However, the SNEC model cannot reproduce the light curves well during the fall from the plateau since the radiation diffusion approach used in SNEC is no longer valid during and after this period, so we were not able to use SNEC to explore the effect of ⁵⁶Ni mixing on the postplateau light curves. Morozova et al. (2015) also found that the light curves generated by SNEC are not sensitive to the degree of ⁵⁶Ni mixing, so we fixed the initial ⁵⁶Ni mixing and mixed ⁵⁶Ni up to 5 M_⊙ in the mass coordinates.

Morozova et al. (2018) pointed out that only the early phase of the light curve is dominated by CSM, so it is possible to adopt a two-step approach to fit the light curves. In the first step, we evaluated the fit for the part of the light curve that is mostly dominated by the hydrogen-rich envelope and varied only M and E. The fitting range was chosen to be between the end of s1 (37.19 days from explosion for SN 2018cuf) and t_PT, where s1 is the initial steeper slope of the light curve. This allowed us to determine the best-fit progenitor mass and explosion energy. In the next step, we fixed the progenitor mass and explosion energy found in step one, and we explored the influence of CSM, varying K and R_ext and fitting the whole light curve through t_PT. This substantially reduced the number of models needed to explore the parameter space, allowing us to search over a finer grid in each parameter.

At each stage, the best-fit model was determined by interpolating the models to the observed epochs in $g,r,i$ filters and minimizing ${\chi }^{2}$ over these filters. For the first stage, the range of parameters considered is $10\ {M}_{\odot }\ \lt$ M $\lt \,30\ {M}_{\odot }$ and $0.1\lt E\lt 1.2$ (in units of ${10}^{51}\$ erg). We obtained a best fit of $M=14.5\ {M}_{\odot }$, which corresponds to an 827 R_⊙ progenitor star from Table 2 of Sukhbold et al. (2016), and $E=5.71\,\times {10}^{50}\$ erg as shown in the upper panel of Figure 13. In the next step, the CSM parameter range was set to $2\lt K\lt 20$ (in units of ${10}^{17}\,{\rm{g}}\,{\mathrm{cm}}^{-1}$) and $827\ {R}_{\odot }\ \lt$ ${R}_{\mathrm{ext}}\ \lt$ $3000\ {R}_{\odot }$, and the fitting range also included the early part of the light curve, i.e., we fit the light curves from the explosion to t_PT. The result is presented in the middle panel of Figure 13, and the best-fitting model is $K=3.1\times {10}^{17}\,{\rm{g}}\,{\mathrm{cm}}^{-1}$ and ${R}_{\mathrm{ext}}=1369\ {R}_{\odot }$. In the bottom panel of Figure 13, we show the light curves of the best-fitting models with and without dense CSM. The progenitor mass (14.5 M_⊙) we got from the SNEC model is in good agreement with what we got from the synthetic nebular spectrum analysis (12–15 M_⊙), and this is a moderate mass for an SN II. It should be noted that we did not fit the model photospheric velocity with the ejecta velocity derived from the [Fe ii] λ5169 line. Both Goldberg et al. (2019) and Goldberg & Bildsten (2020) pointed out that fitting the ejecta velocities inferred from the [Fe ii] λ5169 line can barely break the degeneracy between the explosion properties, so we chose not to fit the ejecta velocity in our SNEC modeling.

Zoom In Zoom Out Reset image size

Previous work (Morozova et al. 2017, 2018; Bostroem et al. 2019) has found that there is a strong degeneracy between the density profile and the external radius of CSM and the total mass of CSM derived from the fits is more robust. If we adopt a progenitor radius of 827 R_⊙ as the inner CSM radius, the total CSM mass of our best-fit model will be 0.07 M_⊙ by integrating Equation (6) over r. If we interpret this wind as that of a typical RSG, adopting a wind speed of 10 km s⁻¹, the mass-loss rate would be 0.06 M_⊙ yr⁻¹ within a timescale of 14 months, much higher than those of the steady winds observed in RSGs (Smith 2014). The possible explanation is that such dense CSM may originate from pre-SN outbursts due to the late-stage nuclear burning in the stellar interior (Quataert & Shiode 2012; Smith & Arnett 2014; Fuller 2017; Ouchi & Maeda 2019; Morozova et al. 2020). Due to the presence of dense CSM around the SN, flash signatures may also be seen in the early spectrum (Yaron et al. 2017; Nakaoka et al. 2018; Rui et al. 2019). However, such a signature is not found for SN 2018cuf, which may imply that the dense CSM is very close to the progenitor, consistent with the small CSM radial extent derived from the SNEC model.

After the shock breakout, the SN emission is dominated by shock cooling and carries useful information about the radius and pre-explosion evolution of the progenitor star system. Sapir & Waxman (2017) updated the model presented by Rabinak & Waxman (2011) and found that the photospheric temperature and bolometric luminosity during the early phase for an SN with an RSG progenitor (convective envelope with n = 3/2) can be written as

$$\begin{eqnarray}&&T(t)={T}_{1}\ * \ {t}_{{\rm{d}}}^{-0.45}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&L(t)={L}_{1}\ \exp \left[-{\left(\frac{1.67* {t}_{{\rm{d}}}}{{t}_{\mathrm{tr}}}\right)}^{0.8}\right]\ * \ {t}_{{\rm{d}}}^{-0.17},\end{eqnarray} \tag{ 8 }$$

where T₁ and L₁ are the temperature and the luminosity ∼1 day after the explosion, respectively, t_d is the time from explosion, and t_tr is the time when the envelope is starting to become transparent. We applied this model to SN 2018cuf, which was discovered well before t_tr, using the code developed by Hosseinzadeh et al. (2018) and Hosseinzadeh (2020). This Markov Chain Monte Carlo (MCMC) routine was adopted to give the posterior probability distributions of T₁, L₁, t_tr, and t₀ simultaneously, where t₀ is the explosion epoch. This analytical model is only valid for T < 0.7 eV, and we have checked that the final fitting results satisfy this condition. The MCMC converges to an explosion epoch of MJD 58,287.8 ± 0.2 (or JD 2,458,288.3 ± 0.2), which is about 3 days earlier than our last nondetection (JD 2,458,291.74). An explosion epoch earlier than our first nondetection is possible as the SN may be below our detection limits shortly after explosion. However, in order to fit the U- and V-band light curves, we required t_tr = 10,000 days, which is unphysically late. For this reason we did not attempt to derive progenitor or explosion parameters using this method. The inability of this method to fit the blue part of the light curve has been noted by several authors (Arcavi et al. 2017; Hosseinzadeh et al. 2018). One possible reason for the fitting failure could be that there is a CSM–ejecta interaction around the progenitor, which is supported by the light-curve modeling as we discussed in the last subsection. In addition, the effect of UV-band line blanketing is underestimated in the model spectrum, so that assuming blackbody radiation cannot well reproduce the light curves in UV bands.