Weak Mass Loss from the Red Supergiant Progenitor of the Type II SN 2021yja

SN 2021yja was discovered by ATLAS on 2021-09-08.55 (all dates are in UTC) at c = 15.334 ± 0.007 mag and was not detected on 2021-09-06.48 to a 3σ limit of o > 19.22 mag (Smith et al. 2021; Tonry et al. 2021). Its J2000 coordinates, as measured by Gaia Photometric Science Alerts, are α = 03^h24^m21180, $\delta =-21^\circ 33^{\prime} 56\buildrel{\prime\prime}\over{.} 20$ (Hodgkin et al. 2021), 99'' southwest of the center of its host galaxy, NGC 1325. ⁴⁸ It was recovered in earlier unfiltered (357–871 nm) imaging by the DLT40 Survey on 2021-09-08.29 at 16.431 ± 0.028 mag and was not detected on 2021-09-07.28 to a 3σ limit of 18.671 mag. Kilpatrick (2021) also recovered it at g = 20.4 ± 0.2 mag and r = i = 20.7 ± 0.2 mag in publicly available imaging taken on 2021-09-07.63 with the Multicolor Simultaneous Camera for Studying Atmospheres of Transiting Exoplanets 3 (MuSCAT3) on Las Cumbres Observatory's 2 m Faulkes Telescope North (Narita et al. 2020), taken as part of the Faulkes Telescope Project's education and public outreach program. ⁴⁹ The very low luminosities implied by these magnitudes meant that it was not clear in real time whether these MuSCAT3 images were taken before or after explosion (see Section 3.2).

SN 2021yja was initially classified as a young SN II by Pellegrino et al. (2021) using a spectrum taken by the Global Supernova Project (GSP) on 2021-09-09.51, 15 hours after the ATLAS discovery report. Their classification was based on spectroscopic comparisons with other SNe II using Superfit (Howell et al. 2005) and the Supernova Identification code (SNID; Blondin & Tonry 2007). Shortly afterward, it was reclassified as a young stripped-envelope SN by Williamson et al. (2021) using a second spectrum taken by the GSP on 2021-09-09.60. They argued that the spectrum also resembled other young SN Ic spectra, and that the putative high-velocity Hα line could have been misidentified silicon 635.5 nm. Several days later, after stronger hydrogen lines had developed, Deckers et al. (2021a, 2021b) restored the SN II classification using a spectrum taken on 2021-09-14.31 by the advanced extended Public ESO Spectroscopic Survey of Transient Objects (ePESSTO+; Smartt et al. 2015).

We obtained UV, optical, and infrared photometry of SN 2021yja using the Sinistro cameras on Las Cumbres Observatory's network of 1 m telescopes (Brown et al. 2013), ⁵⁰ MuSCAT3 on Las Cumbres Observatory's 2 m Faulkes Telescope North (Narita et al. 2020), an Andor iKON-L 936 BV camera on the 0.7 m telescope at Thacher Observatory (Swift et al. 2022), KeplerCam on the 1.2 m telescope at F. L. Whipple Observatory (Szentgyorgyi et al. 2005), the Lulin Compact Imager on the 1 m telescope at Lulin Observatory, the MMT and Magellan Infrared Spectrograph (MMIRS) on the MMT (McLeod et al. 2012), the Direct Imaging Camera on the Nickel Telescope (Stone & Shields 1990), the Ultraviolet/Optical Telescope on Swift (Roming et al. 2005), and Andor Apogee Alta F47 cameras on the Panchromatic Robotic Optical Monitoring and Polarimetry Telescopes at Cerro Tololo Inter-American Observatory and Meckering Observatory (Reichart et al. 2005) as part of the DLT40 Survey. See Appendix A for details of the photometry reduction. We also include publicly available photometry from the ATLAS forced photometry server (Tonry et al. 2018a; Smith et al. 2020) and the Supernova Observations and Simulations group (Martinez et al. 2021). Figure 1 shows the light curve of SN 2021yja; the data are available in machine-readable format in the online journal.

Zoom In Zoom Out Reset image size

In addition, we obtained high-resolution follow-up imaging with the Gemini South Adaptive Optics Imager (GSAOI) at Cerro Pachón, Chile (McGregor et al. 2004). We observed the site of SN 2021yja on 2021-12-18 in H-band using an on-off pattern for 300 s of total on-source exposure. Following procedures in Kilpatrick et al. (2021), we performed image calibration using dark and flat-field frames obtained in the instrumental configuration.

After classification, we initiated a high-cadence optical spectroscopic follow-up campaign using the FLOYDS instruments on Las Cumbres Observatory's 2 m Faulkes Telescopes North and South (FTN & FTS; Brown et al. 2013), Binospec on the MMT (Fabricant et al. 2019), the Wide Field Spectrograph (WiFeS) on the 2.3 m telescope at Siding Spring Observatory (SSO; Dopita et al. 2007, 2010), the ESO Faint Object Spectrograph and Camera 2 (EFOSC2) on the New Technology Telescope (NTT; Buzzoni et al. 1984); the Keck Cosmic Web Imager (KCWI) on Keck II (Morrissey et al. 2018), the Robert Stobie Spectrograph (RSS) on the South African Large Telescope (SALT; Smith et al. 2006), the Kast Spectrograph on the Shane Telescope (Miller & Stone 1994), the Boller & Chivens Spectrograph (B&C) on the Bok Telescope (Green et al. 1995), the High-Resolution Echelle Spectrometer (HIRES) on Keck I (Vogt et al. 1994), the Goodman High Throughput Spectrograph on the Southern Astrophysical Research Telescope (SOAR; Clemens et al. 2004), and the Dual Imaging Spectrograph (DIS) on the Astrophysical Research Consortium (ARC) 3.5 m telescope (Lupton 2005). Figure 2 compares the earliest spectra of SN 2021yja, taken within 6 days of estimated explosion, to other early spectra of SNe II from Yaron et al. (2017); Andrews et al. (2019); Hiramatsu et al. (2021b); Tartaglia et al. (2021), and Terreran et al. (2022). The remaining optical spectra are plotted in Figure 3. We also obtained near-infrared spectra using the Near-Infrared Echelette Spectrometer (NIRES) on Keck II (Wilson et al. 2004), the Son of ISAAC spectrograph (SOFI) on the NTT (Moorwood et al. 1998), TripleSpec 4.1 on SOAR (Schlawin et al. 2014), and MMIRS on the MMT, which are plotted in Figure 4. All spectra were either observed at the parallactic angle, or an atmospheric dispersion corrector was used. Details of the spectroscopic reductions are given in Appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

NGC 1325 and the site of SN 2021yja were observed by HST on 1997-03-26 with the Wide Field Planetary Camera 2 (WFPC2) in F606W (PI Stiavelli, SNAP-6359). We reduced all c0m frames with the hst123 reduction pipeline and performed final photometry in the c0m frames using dolphot. In addition, we obtained DECam grizY imaging taken between 2014 and 2019 and covering the site of SN 2021yja. We stacked and calibrated all such data using a custom pipeline based on the photpipe imaging and photometry package (Rest et al. 2005). Finally, we downloaded pre-explosion Spitzer/IRAC imaging (Fazio et al. 2004; Werner et al. 2004; Gehrz et al. 2007) from the Spitzer Heritage Archive and combined the basic calibrated data (cbcd) frames following methods described in Kilpatrick et al. (2021). These pre-explosion data are analyzed in Section 3.5.

NGC 1325 has a heliocentric redshift of z = 0.005307 ±0.000005 (Springob et al. 2005). However, we observe a weak, narrow Hα emission line in our Keck HIRES spectrum at a redshift of z = 0.00568, corresponding to a line-of-sight velocity of cΔz = + 112 km s⁻¹ with respect to the galaxy core. This velocity is reasonable given the position of the SN and a typical galaxy rotation curve (Falcón-Barroso et al. 2017; Guidi et al. 2018), so we adopt this as the SN redshift.

The equivalent widths of Na i D absorption lines have been shown to correlate with extinction due to interstellar dust (Richmond et al. 1994; Munari & Zwitter 1997). In order to measure the extinction in the direction of SN 2021yja, we examine the Na i D features in our Keck HIRES spectrum (pixel scale Δx ≈ 0.0026 nm in the regions of interest), shown in Figure 5. The spectrum reveals four distinct absorption systems at z = 0.00002 (Milky Way), 0.00559, 0.00566, and 0.00568.

Zoom In Zoom Out Reset image size

We measure the equivalent widths of each of these eight lines by fitting and integrating a Gaussian profile, resulting in the values shown in Table 1. (For the systems at z = 0.00566 and 0.00568, we simultaneously fit two Gaussian profiles.) We then convert these equivalent widths to E(B − V) using Eq. (9) of Poznanski et al. (2012) and applying the renormalization factor of 0.86 from Schlafly et al. (2010). The Milky Way extinction value, $E{(B-V)}_{\mathrm{MW}}={0.0152}_{-0.0026}^{+0.0031}$ mag, is consistent with the value from Schlafly & Finkbeiner (2011), E(B − V)_MW = 0.0191 mag. We adopt the latter. The host galaxy extinction value (from the sum of the six lines) is $E{(B-V)}_{\mathrm{host}}={0.085}_{-0.014}^{+0.017}$ mag.

This is a surprisingly high value given the already very blue observed colors of SN 2021yja (see Section 3.1). To confirm this measurement, we also measure the equivalent width of the diffuse interstellar band (DIB) at 578 nm. From Equation (6) of Phillips et al. (2013), which can be rewritten

$$\begin{eqnarray}&&{W}_{\lambda }=({0.059}_{-0.024}^{+0.039}\,\mathrm{nm}\,{\mathrm{mag}}^{-1})E(B-V)\end{eqnarray} \tag{ 1 }$$

using the extinction law of Schlafly & Finkbeiner (2011), we would expect an equivalent width of ${W}_{\lambda }={0.0050}_{-0.0020}^{+0.0033}$ nm. Due to the lower signal-to-noise ratio of this feature, we fit a Gaussian with a fixed center (578 nm at z = 0.00568) and width (FWHM = 0.222 nm; Tuairisg et al. 2000). Again integrating the best-fit Gaussian, we get an equivalent width consistent with our expectation, W_λ = 0.0076 nm. Given the more robust detection of the Na i D lines, we adopt the extinction derived above, for a total of E(B − V) = 0.104 mag, for which we correct using the extinction law of Fitzpatrick (1999). Note that this only accounts for interstellar extinction in the host galaxy, and not for any extinction due to CSM.

Distance estimates for NGC 1325 based on the method of Tully & Fisher (1977) listed on the NASA/IPAC Extragalactic Database range from 17.7 to 26.1 Mpc (μ = 31.24 − 32.09 mag). For comparison, the redshift-dependent distance of NGC 1325 is 23.6 ± 4.5 Mpc, assuming the cosmological parameters of the Planck Collaboration et al. (2020) and an uncertainty of cz ≈ 300 km s⁻¹ due to the galaxy's peculiar velocity.

As an additional check on the SN distance, we apply the expanding photosphere method (EPM; Kirshner & Kwan 1974), which is a geometrical technique that can independently constrain the distance of an individual Type II SN. Assuming that the SN photosphere is expanding freely and spherically shortly after the explosion, we can derive 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}_{{phot}}}\right)+{t}_{0}\end{eqnarray} \tag{ 2 }$$

where D is the distance, v_phot and θ are the velocity and angular radius of the photosphere, respectively, and t₀ is the explosion epoch. We estimate the photospheric velocity by measuring the minimum of the P Cygni profile of the Fe ii line at 516.9 nm, which becomes detectable ∼ 19 days after explosion. These measurements are listed in Table 2. In practice, the SN photosphere radiates as a dilute blackbody, so a dilution factor has to be involved when determining θ (Eastman et al. 1996; Hamuy et al. 2001; Dessart & Hillier 2005). We combine the multiband photometry to simultaneously derive the angular size (θ) and color temperature (T_c) by minimizing the equation

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

where ξ and b_ν are the dilution factor and the synthetic magnitude, respectively, both of which 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 subset, i.e., {B, V}, {B, V, I} and {V, I}. We use the dilution factors calculated by Jones et al. (2009) with the atmosphere model from Dessart & Hillier (2005). The best-fitting parameters D and t₀ in Equation (2) are estimated with a Markov Chain Monte Carlo (MCMC) method, and the priors of parameters are assumed to be uniform. The upper bound on the prior of t₀ is set to be MJD 59464.63, the first detection of the SN. Jones et al. (2009) showed that there is a clear departure from linearity between θ/v and t ≳ 40 days after explosion. Therefore, only data before this phase are used in our calculations (see Table 2). For the three filter subsets, we obtain distances of ${24.0}_{-1.4}^{+2.5}$ Mpc, ${22.4}_{-0.5}^{+0.8}$ Mpc, and ${23.5}_{-1.3}^{+2.2}$ Mpc, respectively. The explosion epochs are constrained to be $-{2.2}_{-3.2}^{+1.7}$ d, $-{0.6}_{-1.0}^{+0.6}$ d, and $-{2.1}_{-3.0}^{+1.6}$ d, respectively, relative to the estimate from our shock cooling model (see Section 3.2).

Table 2. Data and Results for the Expanding Photosphere Method

Phase	Fe ii Velocity	B	V	I	BV		BVI		VI
(d)	(Mm s−1)	(mag)	(mag)	(mag)	θ/v (Mpc−1 d)	Temp. (K)	θ/v (Mpc−1 d)	Temp. (K)	θ/v (Mpc−1 d)	Temp. (K)
19.0	8.122 ± 0.389	14.733 ± 0.056	14.549 ± 0.062	14.256 ± 0.045	0.805 ± 0.127	${11533}_{-1500}^{+2103}$	0.847 ± 0.068	${10735}_{-573}^{+640}$	0.866 ± 0.123	${10438}_{-1004}^{+1308}$
22.0	8.091 ± 0.161	14.894 ± 0.061	14.582 ± 0.054	14.201 ± 0.029	0.940 ± 0.100	${9353}_{-989}^{+1260}$	0.990 ± 0.044	${9213}_{-357}^{+395}$	1.015 ± 0.087	${9258}_{-632}^{+759}$
23.0	7.747 ± 0.180	14.894 ± 0.061	14.548 ± 0.063	14.192 ± 0.028	1.032 ± 0.113	${8907}_{-942}^{+1228}$	1.040 ± 0.047	${9212}_{-356}^{+394}$	1.028 ± 0.103	${9569}_{-758}^{+922}$
24.0	7.599 ± 0.161	14.901 ± 0.056	14.518 ± 0.058	14.178 ± 0.027	1.105 ± 0.104	${8465}_{-788}^{+989}$	1.075 ± 0.045	${9138}_{-333}^{+362}$	1.032 ± 0.100	${9767}_{-744}^{+898}$
25.9	7.190 ± 0.139	15.012 ± 0.090	14.605 ± 0.067	14.213 ± 0.040	1.131 ± 0.136	${8349}_{-1035}^{+1412}$	1.153 ± 0.063	${8757}_{-447}^{+500}$	1.146 ± 0.124	${9170}_{-785}^{+975}$
26.2	7.264 ± 0.142	15.022 ± 0.076	14.609 ± 0.066	14.210 ± 0.040	1.131 ± 0.123	${8211}_{-906}^{+1192}$	1.154 ± 0.059	${8646}_{-402}^{+445}$	1.145 ± 0.121	${9093}_{-761}^{+942}$
28.1	7.220 ± 0.138	15.049 ± 0.076	14.614 ± 0.057	14.223 ± 0.037	1.159 ± 0.111	${7968}_{-806}^{+1039}$	1.155 ± 0.056	${8644}_{-383}^{+424}$	1.136 ± 0.109	${9166}_{-696}^{+844}$
29.7	5.986 ± 0.177	15.148 ± 0.073	14.643 ± 0.054	14.170 ± 0.012	1.456 ± 0.125	${7333}_{-655}^{+815}$	1.505 ± 0.059	${7965}_{-263}^{+324}$	1.535 ± 0.106	${8368}_{-485}^{+638}$
30.1	6.719 ± 0.115	15.166 ± 0.065	14.642 ± 0.054	14.169 ± 0.013	1.313 ± 0.101	${7189}_{-593}^{+709}$	1.353 ± 0.037	${7861}_{-232}^{+249}$	1.368 ± 0.090	${8370}_{-488}^{+649}$
31.2	6.312 ± 0.111	15.240 ± 0.080	14.658 ± 0.058	14.181 ± 0.018	1.439 ± 0.122	${6791}_{-613}^{+756}$	1.449 ± 0.048	${7728}_{-275}^{+322}$	1.452 ± 0.104	${8340}_{-521}^{+699}$
32.1	6.041 ± 0.093	15.202 ± 0.076	14.625 ± 0.059	14.193 ± 0.026	1.524 ± 0.126	${6813}_{-596}^{+734}$	1.485 ± 0.052	${7924}_{-292}^{+317}$	1.447 ± 0.114	${8723}_{-592}^{+710}$
33.9	6.093 ± 0.096	15.317 ± 0.070	14.696 ± 0.056	14.131 ± 0.033	1.508 ± 0.114	${6499}_{-508}^{+606}$	1.613 ± 0.054	${7131}_{-233}^{+252}$	1.670 ± 0.101	${7578}_{-440}^{+511}$
37.0	5.659 ± 0.089	15.401 ± 0.061	14.684 ± 0.051	14.162 ± 0.029	1.734 ± 0.116	${5938}_{-384}^{+444}$	1.729 ± 0.052	${7045}_{-202}^{+214}$	1.713 ± 0.102	${7908}_{-441}^{+498}$
39.8	5.553 ± 0.075	15.516 ± 0.067	14.757 ± 0.054	14.195 ± 0.028	1.749 ± 0.125	${5736}_{-391}^{+448}$	1.766 ± 0.049	${6809}_{-197}^{+212}$	1.776 ± 0.098	${7598}_{-420}^{+475}$

A machine-readable version of the table is available.

Download table as:  DataTypeset image

Seeing that these distances are consistent with the Tully-Fisher distance above, and given the large uncertainties from this method, we adopt the most recent Tully-Fisher estimate (calculated using the I band), ${23.4}_{-4.4}^{+5.4}$ Mpc (μ = 31.85 ± 0.45 mag; Tully et al. 2016). However, we keep in mind that all of our luminosity-dependent measurements suffer from a large systematic uncertainty in the distance.

For the purposes of constructing a bolometric light curve, we restrict ourselves to epochs where we have Swift observations, as UV photometry is critical for constraining the SED of young SNe. We supplement the Swift photometry with ground-based griz observations taken within 0.5 days of the Swift observation window. We fit a blackbody spectrum to the observed SED at each epoch using an MCMC routine implemented in the Light Curve Fitting package (Hosseinzadeh & Gomez 2020). This gives the temperature and radius evolution in the bottom panel of Figure 6. We then integrate the best-fit blackbody from the blue edge of the U band to the red edge of the I band to produce a pseudobolometric light curve, such that it will be comparable to pseudobolometric light curves in the literature. The top panel of Figure 6 compares the pseudobolometric light curve of SN 2021yja to 39 other SNe II from Valenti et al. (2016). SN 2021yja is among the most luminous.

Zoom In Zoom Out Reset image size

A comparison with the light-curve model grid from Hiramatsu et al. (2021a) suggests that a hydrogen-rich envelope mass ≳ 5 M_☉ is required to reproduce the unusually long plateau duration ( ≈ 140 days). A higher explosion energy ( ≳ 2 × 10⁵¹ erg) is also likely required to reproduce the plateau luminosity of SN 2021yja, which is higher than the models, suggesting an even more massive hydrogen-rich envelope according to light-curve scaling relations (e.g., Popov 1993; Kasen & Woosley 2009; Goldberg et al. 2019).

As suggested by Rubin & Gal-Yam (2017) and Kozyreva et al. (2020), the temperature evolution prior to the recombination phase can serve as a solid diagnostic of the progenitor radius. The temperature declines as a power law during this period and changes to a shallower decline after the onset of recombination (Nakar & Sari 2010; Shussman et al. 2016; Sapir & Waxman 2017; Faran et al. 2019), making this phase relatively easy to distinguish. In the bottom panel of Figure 6, we compare the temperature evolution of SN 2021yja to temperature models from Kozyreva et al. (2020) for various progenitor radii. We find that SN 2021yja lies between the 800 R_☉ and 1000 R_☉ models. Although there is no larger/hotter model, this comparison demonstrates that the progenitor of SN 2021yja is at least consistent with typical RSG radii (100 − 1500 R_☉; Levesque 2017), in contrast to our findings in Section 3.2.

We also downloaded UV light curves of all SNe II in the Swift Optical/Ultraviolet Supernova Archive (Brown et al. 2014). We correct these for Milky Way extinction only and plot them in Figure 7. SN 2021yja is the second most UV-luminous SN II observed by Swift, after only SN 2020pni (Terreran et al.2022).

Zoom In Zoom Out Reset image size

Lastly, in Figure 8, we compare the optical colors of SN 2021yja to a sample of other SNe II (Hamuy et al. 2001; Leonard et al. 2002; Elmhamdi et al. 2003; Bose et al. 2013; Munari et al. 2013; Anderson et al. 2014; Brown et al. 2014; Dall'Ora et al. 2014; Faran et al. 2014; Galbany et al. 2016; Rubin et al. 2016; Hosseinzadeh et al. 2018; de Jaeger et al. 2019; Szalai et al. 2019; Andrews et al. 2019; Bostroem et al. 2020; Dong et al. 2021; Terreran et al. 2022). The evolution is fairly typical, starting near the minimum colors for a hot blackbody, then evolving redward, with an inflection point around 15 days after explosion. In the terminology of Anderson et al. (2014), this is the transition time t_tran when the light-curve slope changes from s₁ (the initial decline) to s₂ (the "plateau"; see their Figure 1). Notably, we observe this transition in the colors even though the bolometric light curve in Figure 6 does not have a significant inflection point. We discuss this further in Section 4.1.

Zoom In Zoom Out Reset image size

In the absence of circumstellar interaction, the early light curves of CCSNe are dominated by shock cooling emission, which depends upon the parameters of the explosion and the progenitor star, in particular the progenitor radius. We fit the model of Sapir & Waxman (2017) to the early light curve of SN 2021yja using the MCMC routine implemented in the Light Curve Fitting package (Hosseinzadeh & Gomez 2020). Because of the large uncertainties in extinction and distance outlined above, we allowed these to be free parameters in our model, with Gaussian priors according to their known values. We also included an intrinsic scatter term, σ, such that the effective uncertainty on each point is increased by a factor of $\sqrt{1+{\sigma }^{2}}$, with a half-Gaussian prior. This allows us to capture additional uncertainties, e.g., from photometric calibration, deviations from a blackbody spectrum, or small variations around the model, in our reported parameter uncertainties. The prior shapes and limits are given in Table 3. We ran 100 walkers for 1000 steps to reach convergence, assessed via visual inspection of the chain history, and for an additional 1000 steps to sample the posterior. Figure 9 shows the best-fit light curve, as well as the posterior distributions of and correlations between each parameter.

Zoom In Zoom Out Reset image size

Table 3. Shock Cooling Parameters

Parameter	Variable	Prior Shape	Prior	Parameters a	Best-fit Value b	Units
Shock velocity	vs*	Uniform	0	3	1.0 ± 0.2	108.5 cm s−1
Envelope mass c	Menv	Uniform	0	10	0.6 ± 0.1	M☉
Ejecta mass × numerical factor d	fρ M	Uniform	3	100	60 ± 30	M☉
Progenitor radius	R	Uniform	0	100	14 ± 2	1013 cm
Distance	dL	Gaussian	23.4	4.9	25 ± 5	Mpc
Extinction	E(B – V)	Gaussian	0.104	0.0155	${0.105}_{-0.009}^{+0.010}$	mag
Explosion time	t0	Uniform	59460	59464.6	59464.40 ± 0.06	MJD
Intrinsic scatter	σ	Half-Gaussian	0	1	4.8 ± 0.2	Dimensionless

Notes.

^a The "Prior Parameters" column lists the minimum and maximum for a uniform distribution, and the mean and standard deviation for a Gaussian distribution. ^b The "Best-fit Value" column is determined from the 16th, 50th, and 84th percentiles of the posterior distribution, i.e., median ± 1σ. ^c Sapir & Waxman (2017) define the envelope of an RSG as the region where $\delta \equiv \tfrac{R-r}{R}\ll 1$, where the density follows ${\rho }_{0}(\delta )=\tfrac{3{f}_{\rho }M}{4\pi {R}^{3}}{\delta }^{n}$ (we adopt $n=\tfrac{3}{2}$ for convective envelopes). ^d The ejecta mass and density profile do not have a strong effect on the light curve. Therefore this parameter is essentially unconstrained.

Download table as:  ASCIITypeset image

Qualitatively, the model fits SN 2021yja reasonably well across all bands. Unlike for several recent SNe II that have a significant peak in their early light curve (see, e.g., Hosseinzadeh et al. 2018), the shock cooling model is able to reproduce the smooth rise to plateau over ∼ 10 days in the reddest bands. The model is nominally valid until MJD 59496.3 (31.7 days after explosion), at which point the temperature drops below 0.7 eV and recombination breaks the assumption of constant opacity (Sapir & Waxman 2017). This phase is roughly consistent with the recombination limit we estimated in Figure 6. However, there are clear hydrogen features in the spectra starting around 10 days, implying that recombination is occurring much earlier than this. For this reason, we do not fit our light curve after maximum light in the reddest bands (i.e., we only fit the data shown in Figure 9). This inconsistency in recombination time is one clue that the physics behind the shock cooling model may not completely describe the data, despite the appearance of a good fit.

Putting this concern aside for the moment, the model is most sensitive to the progenitor radius. The best-fit progenitor radius is ${2030}_{-230}^{+260}\ {R}_{\odot }$, which is significantly larger than typical RSG radii (100 − 1500 R_☉; Levesque 2017), and larger than the radius we infer in Section 3.1. Inspecting the correlations in Figure 9, we see that the radius estimate is strongly (positively) correlated with the extinction estimate. In other words, if we decrease our adopted extinction by 1 − 2σ, we would estimate a more reasonable radius. However, we performed a second fit with a uniform prior in E(B − V) and derived a best-fit extinction value similar to that in Section 2.5. (As Sapir & Waxman 2017 assume a blackbody SED in their model, this is equivalent to determining the extinction required to redden a blackbody into the observed SED.) Alternatively, as we discuss further in Section 4.1, the extreme UV flux may be the effect of circumstellar interaction on our light curve, which is not included in the shock cooling model. If we could model this interaction component separately, a smaller radius might be required to generate the shock cooling component. We can treat our best-fit radius as a robust upper limit on the true progenitor radius, though admittedly it is not very constraining.

We also find a best-fit explosion (shock-breakout) time of MJD 59464.40 ± 0.06, a mere 5.4 ± 1.4 rest-frame hours before our serendipitous detection with MuSCAT3. One possibility, first suggested by Kilpatrick (2021), is that this first detection is of the progenitor prior to shock breakout. Indeed the absolute magnitude at this phase (M_r = − 11.5 mag, M_i = − 11.3 mag) is similar to the "precursor emission" of SN 2020tlf observed by Jacobson-Galán et al. (2022). To investigate this possibility, we performed the light-curve fit again without this first detection. The best-fit explosion time is still before the MuSCAT3 detection, and the model has a slightly larger intrinsic scatter (σ ≈ 5.8). In the absence of extraordinary proof to the contrary, we conclude that this detection was after explosion and adopt the explosion time above as phase = 0 throughout our figures.

Near the end of our observations presented here, as of around 138 rest-frame days after explosion, the light curve of SN 2021yja has settled on the radioactive-decay tail, where the ⁵⁶Co → ⁵⁶Fe decay provides most of the input power. We confirm that the bolometric decline rate over the last four epochs of photometry from Las Cumbres Observatory is 0.0096 mag day⁻¹, very close to the ⁵⁶Co → ⁵⁶Fe decay rate (0.0098 mag day⁻¹; Colgate & McKee 1969). During this phase, we can measure the mass of radioactive ⁵⁶Ni produced in the explosion by scaling the luminosity to that of SN 1987A, which had an independent ⁵⁶Ni measurement (M_Ni,87A=0.075 M_☉; Spiro et al. 2014; Valenti et al. 2016).

We calculate the bolometric luminosity for these four epochs by directly integrating the SED over the U to i filters. This is different than the blackbody-fitting method we used during the photospheric phase, because at these late times the SED is no longer described by a Planck (1906) function (e.g., Martinez et al.2022). We calculate a bolometric light curve of SN 1987A with the same method using photometry from Catchpole et al. (1988) via the Open Supernova Catalog (Guillochon et al. 2017), interpolate to the four epochs of SN 2021yja photometry, and take the average ratio. We find that SN 2021yja produced ${M}_{\mathrm{Ni}}={0.141}_{-0.049}^{+0.074}\ {M}_{\odot }$, where the uncertainty is dominated by our distance uncertainty. This places SN 2021yja at the 97th percentile for nickel mass among SNe II (Anderson 2019), although we caution that this measurement is only ≈ 1.6σ above the nickel mass of SN 1987A. Another possibility is that our measurement is inflated by extra luminosity from ongoing circumstellar interaction, although the decline slope we measure above suggests that ⁵⁶Co decay is the dominant power source. An unusually large nickel mass may suggest a massive progenitor, according to the correlations of Eldridge et al. (2019). We discuss this possibility further in Section 4.2.

The spectroscopic evolution of SN 2021yja ≳ 6 days after explosion is typical for a SN II: P Cygni lines of hydrogen and helium superimposed on a cooling blackbody. Metal lines, including Ba ii and Fe ii begin to appear several weeks after explosion, blanketing the blue side of the spectrum. Its infrared spectra (Figure 4) place it in the "weak" class of Davis et al. (2019), with weak helium absorption at 1.083 μm in our > 100 d spectra, a high-velocity component of the same line, and more pronounced Sr ii absorption. As such, SN 2021yja is consistent with the correlation between slow plateau decline rates and weak infrared helium absorption (Davis et al. 2019). No carbon monoxide emission is detected at 100 days after explosion.

However, our earliest optical spectra, taken around 2 days after explosion, are somewhat unusual in the sense that they do not show the narrow high-ionization lines seen in many SNe II this early. These lines indicate short-lived circumstellar interaction, which ends when the nearby circumstellar material (CSM) is swept up by the expanding ejecta (Gal-Yam et al. 2014; Smith et al. 2015; Khazov et al. 2016; Yaron et al. 2017; Bullivant et al. 2018; Hosseinzadeh et al. 2018; Soumagnac et al. 2020; Bruch et al. 2021; Hiramatsu et al. 2021b; Terreran et al. 2022). Figure 2 compares the early spectra of SN 2021yja to other SNe II with well-constrained explosion dates and high-quality spectra taken within 6 days of explosion. Some show narrow He ii, C iii, N iii, and other lines: SN 2013fs (Yaron et al. 2017; Bullivant et al. 2018), SN 2017ahn (Tartaglia et al. 2021), SN 2018zd (Hiramatsu et al. 2021b), and SN 2020pni (Terreran et al. 2022). SN 2017gmr (Andrews et al. 2019) and SN 2021yja do not show these lines, but instead show a unique "ledge-shaped" feature around 450–480 nm, which we discuss below.

Figure 10 shows a detailed view of the ledge-shaped feature in the earliest spectrum of SN 2021yja, compared to Hα and He 587.6 nm in the same spectrum. There are two lines of thought about the identity of this feature in the literature, both related to circumstellar interaction. Bullivant et al. (2018, their Figure 20) and Andrews et al. (2019, their Figure 18) interpret it as very broad, blueshifted He ii 468.6 nm (although there is also a narrow component of this line present) produced in the outermost layers of the SN ejecta. However, the line profile is somewhat more symmetrical in these SNe than in SN 2021yja. Soumagnac et al. (2020, their Figure 7) and Bruch et al. (2021, their Figure 5) interpret this feature as a blend of various flash-ionized lines from the CSM: N v 460.4 and 462.0 nm, N ii 463.1 and 464.3 nm, C iv 465.8 nm, and He ii 468.6 nm.

Zoom In Zoom Out Reset image size

Neither of these explanations is fully satisfactory. The former (single, broad line) is difficult to reconcile with the irregular profile of the feature, when other features have typical P Cygni profiles. In addition, the implied velocity is nearly 30,000 km s⁻¹ (as measured from the He ii rest wavelength to the blue edge of the feature), which is much higher than any other feature in the spectrum. The latter explanation (blend of several narrow emission lines) is inconsistent with the breadth and strength of the feature. We suspect that the feature is indeed a blend of several lines, possibly including He ii, C iii, and N iii, although these must be either much broader or much more numerous than the features suggested by previous work. We discuss one such possibility in Section 4.1.

Figure 11 (left) shows our high-resolution GSAOI image of SN 2021yja, aligned to the pre-explosion WFPC2 F606W image using three common astrometric sources in both frames, resulting in an rms astrometric uncertainty of 0.08'' (0.8 WFPC2 pixels) between both frames. In addition, we aligned the GSAOI image to DECam and Spitzer/IRAC imaging, using 4–6 common astrometric standards in both frames and resulting in an rms uncertainty of 0.07–0.12'' (Figure 11, right).

Zoom In Zoom Out Reset image size

The position of SN 2021yja does not align with any counterpart in the HST imaging. The closest source is 0.88'' (11σ) away from this position and detected at m_F606W = 24.05 ± 0.06 mag (previously discussed in Kilpatrick 2021). The same is true in each of the other pre-explosion images as shown in Figure 11, with no evidence for any significant, point-like emission within at least 3σ of the site of SN 2021yja.

Motivated by the depth of our imaging set, we place 3σ limits on the presence of a counterpart in each of our pre-explosion images. We estimate the standard deviation of the background variation within 2'' of SN 2021yja and use this value to estimate the maximum total flux that would not be flagged as a 3σ detection within one PSF width at this location. The resulting limits are given in AB mag in Table 4, which we consider the 3σ limits on the presence of any counterpart in each image.

To place these limits in context, we assume a blackbody spectral energy distribution with a given luminosity and show the region of the Hertzsprung-Russell diagram ruled out by our limits assuming a source at 23.4 Mpc and the extinction we derive in Section 2.5. We compare these limits to the expected luminosities and temperatures of single-star models derived from stars in MESA Isochrones & Stellar Tracks (MIST; Choi et al. 2016). We can rule out all stars with initial masses >9 M_⊙ as shown in Figure 12. This analysis is dominated by the WFPC2 limit, which is significantly more constraining for the MIST models than DECam and Spitzer.

Zoom In Zoom Out Reset image size

Thus while we rely on the WFPC2 image from 26 March 1997, 24.5 yr before explosion, for our best constraints on the nature of any progenitor star of SN 2021yja, we can place meaningful constraints on that star at later times when the DECam and Spitzer imaging were obtained. These limits are informative in the scenario where that star was highly variable and had a lower flux in F606W compared to its average flux (i.e., similar to variability in α Ori and the recent detection of a dimming RSG in M51 described in Levesque & Massey 2020; Jencson et al. 2022). We further explore the implications of this scenario below.

Our analysis above gives rise to a somewhat confusing picture with regard to CSM around the progenitor of SN 2021yja. On one hand, we do not see some standard signs of CSM, including pronounced, narrow emission lines in early spectra, and, at first glance, the shock cooling model appears to fit our observed light curve reasonably well. On the other hand, the high UV luminosity, very blue colors, and persistent high temperature (resulting in an overly large progenitor radius from the shock cooling model), as well as the ledge-shaped feature in the early spectra, do suggest some level of early circumstellar interaction. HST UV-optical spectra from Vasylyev et al. (2022) show no sign of circumstellar interaction, but these were observed much later (7 days after explosion) than our early spectra showing the ledge. In addition, the detection of delayed radio emission (Alsaberi et al. 2021) from SN 2021yja following an initial non-detection (Ryder et al. 2021) can be taken as an indication of interaction with a non-uniform CSM, similar to the conclusion of Bostroem et al. (2019) for the SN II ASASSN-15oz. SN 2021yja was also observed in the very high-energy γ-ray domain (E > 100 GeV) by the High Energy Stereoscopic System (H.E.S.S.) Imaging Atmospheric Cherenkov Telescope Array, which will yield complementary constraints on the CSM density (H.E.S.S. Collaboration 2022, in preparation) following the work presented by Abdalla et al. (2019) and the H.E.S.S. Collaboration (2022).

One way to resolve this apparent contradiction is to infer a small amount of CSM: less than is required to produce strong emission lines, but still enough to affect the broadband light curve, especially in the UV. We turn to the literature to constrain the mass-loss rate from both directions. By modeling a sample of early SN II spectra with narrow emission lines, Boian & Groh (2020) find CSM density parameters of order D ∼ 10¹⁵⁻¹⁶ g cm⁻¹, corresponding to mass loss of $\dot{M}\equiv 4\pi {v}_{\infty }D\,\gtrsim {10}^{-4}\ {M}_{\odot }{\mathrm{yr}}^{-1}$, where v_∞ is the terminal wind speed. Our mass-loss rate must therefore be $\dot{M}\lesssim {10}^{-4}\ {M}_{\odot }{\mathrm{yr}}^{-1}$.

Dessart et al. (2017) propose a scenario in which the interaction lines are both broadened and blueshifted, which can potentially explain the ledge-shaped feature we discuss in Section 3.4. This occurs when an RSG with radius R_⋆ = 501 R_☉ explodes into low-density CSM ($\dot{M}={10}^{-6}\ {M}_{\odot }{\mathrm{yr}}^{-1}$). They term this the r1w1 model. In this scenario, narrow emission lines are not seen or disappear within hours of explosion, at which time we had not yet observed the spectrum of SN 2021yja. When they add an extended atmosphere onto the RSG (r1w1h model; scale height H_ρ = 0.3R_⋆), the spectroscopic evolution is similar, but the light curve rises faster and is much more luminous in the UV (reaching M_UVW2 ≈ − 19.5 mag). Such an extended atmosphere only requires a small energy injection during the last year before explosion (Smith & Arnett 2014; Morozova et al. 2020). Figure 13 compares the early spectra of SN 2021yja to the r1w1 and r1w1h models of Dessart et al. (2017). Notably, the extended atmosphere or r1w1h could also lead to an overestimate of the progenitor radius by a shock cooling model that was not designed to account for it. Therefore, we believe this to be the most likely scenario for the progenitor of SN 2021yja: an RSG with an extended atmosphere exploding into low-density CSM.

Zoom In Zoom Out Reset image size

Because of our skepticism about the inferred progenitor radius (Section 4.1), the best constraints we have on the progenitor are (1) the luminosity limit from direct imaging of the field before explosion and (2) the inferred nickel mass from the light-curve tail. Unfortunately, these are also somewhat in conflict with each other: the very strict luminosity limit suggests a very low-mass progenitor (with some caveats, discussed below), whereas the very large nickel mass suggests a massive progenitor ( ≳ 20 M_☉; Eldridge et al. 2019). The latter places SN 2021yja close to the upper limit for SN II progenitors (Smartt 2015; Davies & Beasor 2018, 2020a, 2020b), implying that its progenitor star ought to have M_V = − 5 to −6 mag.

The limits from WFPC2 suggest that M_V > − 4.4 mag, consistent with a progenitor star of initial mass < 9 M_⊙ (and following the single-star models of Choi et al. 2016). We consider this a credible limit on the initial mass of the progenitor star derived using the same distance and extinction assumptions as SN 2021yja above, but we acknowledge that there is some tension with the nature of the progenitor star as inferred from the SN light curve. Our pre-explosion limit is dominated by the single epoch of WFPC2/F606W imaging from 1997, 24.5 years before explosion, which provides the deepest limits. The tension between our pre-explosion constraints on the progenitor star and those derived from SN 2021yja could be partly eased if the star were extinguished or variable on timescales of decades before core collapse.

A −0.6 to −1.4 mag discrepancy between this expectation and our derived limits could be explained by variability in the progenitor star V-band magnitude confined to a single epoch 24.5 yr prior to core collapse, similar to the variability observed in RSGs such as Betelgeuse and M51-DS1 (Levesque & Massey 2020; Jencson et al. 2022). The peak-to-peak variability in Betelgeuse was ≈1.2 mag such that even a fraction of this change in the SN 2021yja progenitor star would completely obscure a <15 M_⊙ star in our WFPC2/F606W imaging. We note, however, that these are extreme cases, and evidence from RSGs observed in M31 (Soraisam et al. 2018) and M51 (Conroy et al. 2018) suggests that the vast majority of RSGs in the expected luminosity range for the SN 2021yja progenitor star do not exhibit such extreme variability. Indeed, assuming that all RSGs follow similar V- and I-band variability to that observed in M51, those observed in Conroy et al. (2018) suggest that <8% of those stars will undergo peak-to-peak dimming at M_I > − 0.6 mag. However, this likelihood is highly uncertain for RSGs as a whole, and in particular, those that are <25 yr from core collapse, which may be even more variable than the overall population (e.g., Jacobson-Galán et al. 2022).

Alternatively, the CSM discussed above could significantly obscure the progenitor star in F606W if it contained a significant dust mass. Assuming Beasor et al. (2020) mass-loss rates for an RSG in this mass range ( ≈ 10⁻⁶ M_☉ yr) and corresponding wind speeds (20–25 km s⁻¹), we consider the implied effect on circumstellar extinction. Following methods in Kilpatrick et al. (2018) for a uniform wind with carbonaceous dust grains, we infer that SN 2021yja would experience a line-of-sight extinction A_V ≈ 0.1 mag, in which case circumstellar extinction could not explain an anomalously low optical luminosity for the SN 2021ya progenitor star. However, we note that this is strongly dependent on the assumed mass-loss rate and wind composition, which could be significantly different.