The Early Light Curve of SN 2023bee: Constraining Type Ia Supernova Progenitors the Apian Way

Despite their utility as cosmological distance indicators, the physics behind Type Ia supernovae (SNe Ia), the explosions of carbon–oxygen white dwarfs, is still poorly understood (see Hoeflich 2017 for a review). In the single-degenerate scenario, the white dwarf accretes material from a main-sequence or giant star (Whelan & Iben 1973), whereas the double-degenerate scenario involves a binary system of two white dwarfs (Iben & Tutukov 1984; Webbink 1984). In addition, the explosion may be a delayed detonation (Khokhlov 1991) or a double detonation (Nomoto 1982; Livne 1990; Woosley & Kasen 2011), and the latter may or may not be dynamically driven (Shen et al. 2018). These models have various, sometimes overlapping, observational signatures (see Jha et al. 2019 for a recent review).

With the advent of high-cadence time-domain surveys, the early light curves (starting ≲1 day after explosion) of SNe Ia have been of particular interest as probes of their progenitor systems and explosion mechanisms. Kasen (2010) predicted that the collision between SN ejecta and a nondegenerate binary companion could lead to a blue excess in the first few days after explosion for ∼10% of viewing angles. Similar bumps can be produced by interaction with circumstellar material (CSM; Piro & Morozova 2016), unusual distributions of radioactive ⁵⁶Ni (Noebauer et al. 2017), and the double-detonation mechanism itself (Polin et al. 2019). However, when examined in detail, these models make different predictions for the strength, duration, and color of the early excess.

The first well-sampled early SN Ia light curve was of SN 2011fe, showing a smooth parabolic rise in g band over the first several days after explosion (Nugent et al. 2011). However, over the past decade, several SNe Ia have shown unusual behavior in various filter bands immediately following explosion, ranging from a very fast rise to an initial decrease in luminosity. Normal to overluminous (99aa-like) SNe Ia with early blue or UV excesses include SN 2012cg (Marion et al. 2016), SN 2017cbv (Hosseinzadeh et al. 2017), iPTF16abc (Miller et al. 2018), SN 2021aefx (Ashall et al. 2022; Hosseinzadeh et al. 2022; Ni et al. 2023a), and SN 2021hpr (Lim et al. 2023). Underluminous SNe Ia with strong UV excesses include iPTF14atg (Cao et al. 2015) and SN 2019yvq (Miller et al. 2020; Siebert et al. 2020; Burke et al. 2021; Tucker et al. 2021). MUSSES16D4 (Jiang et al. 2017) and SN 2018aoz (Ni et al. 2022, 2023b), on the other hand, show early red excesses. Normal Type Ia SN 2018oh (Dimitriadis et al. 2019a; Li et al. 2019; Shappee et al. 2019) and overluminous (03fg-like) Type Ia SN 2021zny (Dimitriadis et al. 2023) also show early excesses in their high-cadence unfiltered light curves from Kepler 2 and the Transiting Exoplanet Survey Satellite (TESS), respectively, although the color of this excess is not known. An excess can be ruled out in three other SNe Ia with equivalent data from Kepler (Olling et al. 2015). It is not clear which, if any, of these groups are physically related. Sample analyses by Burke et al. (2022b, 2022a) and Deckers et al. (2022) show excesses in ∼10% of SNe Ia with early photometry.

Here we present observations of SN 2023bee, another overluminous (99aa-like) SN Ia discovered within ∼1 day of explosion and observed at very high cadence by the Distance Less Than 40 Mpc (DLT40) Survey (Tartaglia et al. 2018). We describe these observations in Section 2. The early light curve shows an excess in the U band during the first ∼3 days after explosion, which we fit with various models in Section 3. We analyze the early spectra of SN 2023bee in Section 4 and use our radio observations to place constraints on CSM in Section 5. In Section 6, we discuss the implications of our analysis on the progenitor systems of SNe Ia in general and SN 2023bee in particular.

2.1. Discovery and Classification

2.2. Photometry

After discovery, we initiated a high-cadence optical–UV follow-up campaign using the Sinistro cameras on Las Cumbres Observatory's network of 1 m telescopes (Brown et al. 2013) as part of the Global Supernova Project, the Panchromatic Robotic Optical Monitoring and Polarimetry Telescopes (Reichart et al. 2005) as part of the DLT40 Survey, and the Ultraviolet/Optical Telescope on the Neil Gehrels Swift Observatory (Roming et al. 2005). We measured point-spread function photometry on the Las Cumbres images and aperture photometry on the other images using dedicated pipelines based on PyRAF (Science Software Branch at STScI 2012), HEAsoft (NASA HEASARC 2014), and Photutils (Bradley et al. 2022). Las Cumbres UBV photometry is calibrated to Rubin et al. (1974) and Landolt (1983, 1992) standard fields observed on the same nights at the same sites. Las Cumbres gri photometry is calibrated to the Pan-STARRS1 3π Catalog (Chambers et al. 2016). DLT40 photometry is calibrated to the AAVSO Photometric All-Sky Survey (Henden et al. 2009), with unfiltered photometry calibrated to the r band. Swift photometry uses the updated zero-points of Breeveld et al. (2010) with the time-dependent sensitivity corrections from 2020. UBV and Swift magnitudes are in the Vega system, and gri magnitudes are in the AB system. Figure 1 shows the full multiband light curve, which is available in machine-readable form.

Zoom In Zoom Out Reset image size

2.3. Spectroscopy

We obtained optical spectra of SN 2023bee using the Robert Stobie Spectrograph (RSS; Smith et al. 2006) on the Southern African Large Telescope (SALT; Buckley et al. 2006), the Goodman High-Throughput Spectrograph (Clemens et al. 2004) on the Southern Astrophysical Research Telescope (SOAR), FLOYDS on Las Cumbres Observatory's Faulkes Telescope South (FTS; Brown et al. 2013), the Boller & Chivens (B&C) Spectrograph on the Bok 2.3 m Telescope (Green et al. 1995), and Binospec on the MMT (Fabricant et al. 2019). We reduced these data using the Image Reduction and Analysis Facility (National Optical Astronomy Observatories 1999), a custom pipeline based on the PySALT package (Crawford et al. 2010), and the Goodman (Torres et al. 2017), FLOYDS (Valenti et al. 2014), and Binospec (Kansky et al. 2019) pipelines. We also downloaded the classification spectrum of Zhai et al. (2023) from the Transient Name Server and include it in our analysis. The spectra are plotted in Figure 2 and logged in Table 1. They are available for download from this Letter and from the Weizmann Interactive Supernova Data Repository (Yaron & Gal-Yam 2012).

Zoom In Zoom Out Reset image size

Table 1. Log of Spectroscopic Observations, Velocities, and Equivalent Widths

MJD	Telescope	Instrument	Si ii 636 nm		C ii 658 nm		Phase
			Vel. (km s−1)	EW (nm)	Vel. (km s−1)	EW (nm)	(d)
59,976.763	LJT	YFOSC	25,220 ± 160	9.8 ± 0.9	22,960 ± 520	3.6 ± 1.0	−15.4
59,976.885	SALT	RSS	24,520 ± 80	10.7 ± 1.5	21,510 ± 250	3.4 ± 1.6	−15.3
59,976.985	SALT	RSS	23,980 ± 80	13.0 ± 1.4	20,790 ± 260	4.1 ± 1.5	−15.2
59,977.564	FTS	FLOYDS	23,050 ± 90	12.3 ± 1.0	21,440 ± 670	2.7 ± 1.1	−14.6
59,978.115	SOAR	Goodman	21,740 ± 90	13.4 ± 1.1	20,540 ± 600	2.9 ± 1.2	−14.1
59,979.002	SALT	RSS	21,050 ± 60	14.1 ± 1.2	⋯	⋯	−13.2
59,979.513	FTS	FLOYDS	20,200 ± 90	13.1 ± 0.8	⋯	⋯	−12.7
59,980.705	FTS	FLOYDS	18,880 ± 110	10.8 ± 0.8	⋯	⋯	−11.5
59,981.641	FTS	FLOYDS	17,930 ± 170	11.6 ± 0.7	⋯	⋯	−10.6
59,983.512	FTS	FLOYDS	15,470 ± 190	8.6 ± 0.8	⋯	⋯	−8.7
59,984.630	FTS	FLOYDS	13,710 ± 170	7.7 ± 0.7	⋯	⋯	−7.6
59,985.648	FTS	FLOYDS	13,150 ± 110	6.2 ± 0.8	⋯	⋯	−6.6
59,987.492	FTS	FLOYDS	12,490 ± 80	5.7 ± 0.9	⋯	⋯	−4.8
59,990.502	FTS	FLOYDS	12,070 ± 60	6.0 ± 0.8	⋯	⋯	−1.8
59,991.604	FTS	FLOYDS	12,150 ± 50	6.1 ± 0.8	⋯	⋯	−0.7
59,994.208	Bok	B&C	12,210 ± 50	6.8 ± 0.6	⋯	⋯	+1.9
60,002.620	FTS	FLOYDS	12,220 ± 40	7.5 ± 0.6	⋯	⋯	+10.3
60,006.598	FTS	FLOYDS	11,970 ± 40	7.4 ± 0.6	⋯	⋯	+14.2
60,012.205	MMT	Binospec	11,710 ± 30	5.3 ± 1.2	⋯	⋯	+19.8
60,012.475	FTS	FLOYDS	11,500 ± 60	6.5 ± 0.7	⋯	⋯	+20.0
60,017.172	MMT	Binospec	11,950 ± 30	2.8 ± 0.8	⋯	⋯	+24.7
60,017.514	FTS	FLOYDS	11,720 ± 50	2.9 ± 0.5	⋯	⋯	+25.0

Note. The initial spectrum from the Lijiang 2.4 m Telescope (LJT) was reported to the Transient Name Server by Zhai et al. (2023). Velocities are derived from the minima of best-fit Gaussians (see Section 4).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

2.4. Extinction and Distance

2.5. Radio Nondetections

We obtained radio observations of SN 2023bee with the Karl G. Jansky Very Large Array (VLA) on 2023-02-09.27, ∼9 days after our last nondetection. The observation block was 1 hour long, with 38.45 minutes on-source time for SN 2023bee. Observations were taken in C-band (4–8 GHz) in the B-configuration of the VLA (DDT: 23A-382, PI: S. Sarbadhicary). The observations were obtained in wide-band continuum mode, yielding 4 GHz of bandwidth sampled by 32 spectral windows, each 128 MHz wide sampled by 2 MHz-wide channels with four polarization products. We used 3C138 as our flux, delay, and bandpass calibrator, and J0909 + 0121 as our complex gain calibrator.

Data were calibrated with the VLA CASA calibration pipeline (version 6.4.1-12), which iteratively flags corrupted measurements, applies corrections from the online system (e.g., antenna positions), and applies delay, flux density, bandpass and complex gain calibrations (McMullin et al. 2007). We then imaged the calibrated visibility data set with tclean in CASA. We used multiterm, multifrequency synthesis as our deconvolution algorithm (set with deconvolver=mtmfs in tclean), which approximates the full 4–8 GHz wide-band spectral structure of the sky brightness distribution as a Taylor-series expansion about a reference frequency (in our case, 6 GHz) in order to reduce frequency-dependent artifacts during deconvolution (Rau & Cornwell 2011). We set nterms=2, which uses the first two Taylor terms to create images of intensity (Stokes-I) and spectral index. We sampled the synthesized beam with 03 pixels and imaged out to 1024 to deconvolve any outlying bright sources and mitigate their sidelobes at the primary beam center. Gridding was carried out with the W-projection algorithm (gridder=wproject) with 16 w-planes. Images were weighted with the Briggs weighting scheme (weighting=briggs) using a robust value of 0 to balance point-source sensitivity with high angular resolution and low sidelobe contamination between sources. The final image has a spatial resolution of 102 × 085 (or 183 × 153 pc, assuming a distance of 37 Mpc), and an rms of about 4.5 μJy beam⁻¹.

No radio source was detected at the site of SN 2023bee in the cleaned, deconvolved 6 GHz image at the 3σ level. The flux at the SN location is 8.3 μJy beam⁻¹, and the rms noise in a circular region ∼28'' across, centered on the SN, is 4.8 μJy beam⁻¹. We therefore assume a flux density upper limit of 22.7 μJy, which is equal to the flux density at the SN location plus 3σ noise to be consistent with previous work (Chomiuk et al. 2012). At a distance of 37 Mpc, this corresponds to an upper limit on 6 GHz luminosity of 3.7 × 10²⁵ ergs s⁻¹ Hz⁻¹ (Figure 3).

Zoom In Zoom Out Reset image size

In addition to the VLA, we also made use of radio observations from the Australia Telescope Compact Array (ATCA). Leung et al. (2023) obtained 3σ upper limits of 84 μJy at 5.5 GHz, 60 μJy at 9.0 GHz, 87 μJy at 16.7 GHz, and 129 μJy at 21.2 GHz within 5 days of the explosion. We also obtained new observations at 5.5 GHz and 9.0 GHz at 85 days past explosion to look for any late-time emission due to denser winds. These data were processed in the manner described by Bufano et al. (2014), and yielded 3σ upper limits of 60 and 30 μJy at 5.5 and 9.0 GHz, respectively. The 5.5 GHz ATCA limits are also plotted in Figure 3, yielding an upper limit on 6 GHz luminosity of 9 × 10²⁵ ergs s⁻¹ Hz⁻¹.

Figure 4 shows the color curves of SN 2023bee compared to other SNe Ia. Like other overluminous SNe with an early U-band excess, namely SNe 2017cbv and 2021aefx, there is a sharp change in slope in the U − B color, and to a lesser extent B − V, in the first few days after explosion. This indicates an initially very blue color, rapidly reddening until it matches the slope of the normal SN Ia color curve (e.g., that of SN 2011fe). In SN 2023bee, this initial reddening phase lasts ∼3 days, the same phase that the excess is visible in the U-band light curve. This indicates an additional blue component not present in all SNe Ia that dominates the emission during the first few days.

Zoom In Zoom Out Reset image size

To investigate the origin of this blue component, we compare the multiband light curve to several models from the literature. First, we consider the model of Kasen (2010) in which the ejecta collide with and are shocked by a nondegenerate binary companion. In this model, the strength and duration of the excess depend mostly on the binary separation a, or equivalently the companion radius R, where a/R = 2–3 if we assume it is in Roche-lobe overflow. We construct a model consisting of the SiFTO SN Ia template from Conley et al. (2008) plus an analytical shock component as described by Kasen (2010), where the latter is suppressed based on the viewing angle according to the approximation of Brown et al. (2012). (See Olling et al. 2015 for an alternative approximation.) In this model, the luminosity scale of the shock component is proportional to the product $a{({{Mv}}^{7})}^{1/4}f(\theta )$, where M is the ejecta mass, v is the transition velocity between density profiles in the outer ejecta, and f(θ) is the fraction of the flux that reaches the observer as a function of the viewing angle. The temperature, and therefore color, of the shock emission also depends on a, so it can be fit for separately. However, it is not practical to constrain the remaining parameters simultaneously. Therefore we fix ${{Mv}}^{7}={M}_{\mathrm{Ch}}{(10,000\,\mathrm{km}\,{{\rm{s}}}^{-1})}^{7}$, i.e., a Chandrasekhar mass of ejecta, the outer layers of which are moving at 10,000 km s⁻¹. We caution that fixing this parameter, which can vary by orders of magnitude due to the large exponent, might cause us to underestimate uncertainties on the other parameters. Off-axis systems with high ejecta masses or velocities look the same as on-axis systems with low ejecta masses or velocities.

Because our focus is on the shock component and not on the physics powering the light-curve peak, we scale each band of the SiFTO model so that it matches the observed maximum of SN 2023bee without preserving the colors of the original template. We also fit for small time shifts in U and i relative to the other bands in order to better match the observed light curve. We only include the first 30 days of data in the fit (until MJD 60,005), because the SiFTO template does not match our data well after this time, specifically the slope of the U-band decline and the secondary infrared maximum. The SiFTO template does not extend into the ultraviolet, so we cannot include our Swift data in the fit. For efficiency, we also exclude the filtered DLT40 photometry, and we bin the remaining points by ∼0.01 d. Finally, we include an intrinsic scatter term σ that effectively inflates the photometric uncertainties by a factor of $\sqrt{1+{\sigma }^{2}}$. We perform the fit using a Markov Chain Monte Carlo (MCMC) routine implemented in the Light Curve Fitting package (Hosseinzadeh et al. 2023a). Table 2 describes the prior and posterior distributions of the fit parameters, and Figure 5 (left) shows the best-fit model compared to our observed light curve.

Zoom In Zoom Out Reset image size

Table 2. Companion-shocking Model Parameters

Parameter	Variable	Prior Shape	Prior Parameters a		Best-fit Value b	Units
Companion-shocking Model (Kasen 2010)
Explosion time	t0	Uniform	59,975.12	59,975.65	${59975.3}_{-0.1}^{+0.2}$	MJD
Binary separation	a	Uniform	0	1	${0.055}_{-0.009}^{+0.010}$	1013 cm
Viewing angle	θ	Uniform	0	180	0±16	deg
SiFTO Model (Conley et al. 2008)
Time of B maximum	${t}_{\max }$	Uniform	59,991.50	59,993.50	59,992.04±0.03	MJD
Stretch	s	Log-uniform	0.5	2	0.988±0.004	dimensionless
Time shift in U	Δ tU	Gaussian	0	1	+0.09±0.05	d
Time shift in i	Δ ti	Gaussian	0	1	+0.29±0.05	d
Combined Model
Intrinsic scatter	σ	Half-Gaussian	0	1	11.0±0.3	dimensionless

Notes.

^a The "Prior Parameters" columns list 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σ, with one exception. The posterior for the viewing angle is approximately half-Gaussian peaking at zero, so we use the 0th and 68th percentiles.

Download table as:  ASCIITypeset image

The model provides an excellent fit to the excess in the U band, and it reveals weaker excesses in all other filters compared to the SiFTO template. The best-fit binary separation implies a companion radius of a few solar radii, i.e., a main-sequence star. The posterior for the viewing angle peaks at zero, with a 68th percentile at 16°, i.e., $\theta$ = 0° ± 16°. Note that the posterior for the binary separation covaries weakly with the posterior for the viewing angle because a larger companion off-axis has the same effect as a smaller companion on-axis. As shown in the bottom left panel of Figure 5, the shock component alone overpredicts the Swift photometry by ∼2 mag in UVM2 and UVW2. This effect was also seen in SNe 2017cbv and 2021aefx and can be attributed to the assumption that the shock component radiates as a blackbody (Hosseinzadeh et al. 2017, 2022). In reality, the shock radiation may be subject to line blanketing when it is reprocessed by the SN ejecta, suppressing the UV flux. Unlike in previous events, the shock model is a reasonable fit to the UVW1 (longest-wavelength UV) photometry, which also shows a change in slope between the second and third observations.

To confirm the excesses at other wavelengths, we repeat our fitting procedure after excluding the U-band data. We find that the posteriors do not change significantly. While there is no requirement that the same mechanism produces the excesses in all bands, a complete explanation requires an effect at all optical to near-UV wavelengths, with a larger effect in the bluer bands. As the early spectra do not display unusual features (see Figure 2), the simplest explanation is an additional hot thermal component to the emission, possibly reprocessed through the SN ejecta. We discuss this further in Section 6.

We also consider the double-detonation models of Polin et al. (2019), in which a shell of surface helium detonates before igniting the core of the white dwarf. Some of these models show an initial blue bump in their light curves. We calculated the U − B and B − V colors of the models and compared them to our observations, applying time offsets from 0 to 3 days at 0.5 day intervals and calculating the reduced χ² (${\bar{\chi }}^{2}\equiv {\chi }^{2}/N$, where N is the number of observed points that overlap the model in time) for each model in the grid. The best-matched model is of a 0.9 M_☉ white dwarf edge-lit by a 0.09 M_☉ helium shell, offset by 0.5 days, with ${\bar{\chi }}^{2}=3342$ and N = 32. This model is shown in the right panel of Figure 5. While the model does reasonably reproduce the early color curve of SN 2023bee, both the U-band and B-band excesses are much stronger than observed. In addition, the peak of the light curve is a very poor match. This is expected: SN 2023bee is overluminous, whereas the double-detonation models of Polin et al. (2019) are underluminous (though see Shen et al. 2021 for overluminous double detonations).

Another model for early light-curve excesses comes from shells or clumps of radioactive ⁵⁶Ni near the surface of the ejecta, providing an extra source of photons that can quickly escape to the observer. Magee & Maguire (2020) produce a set of such models designed to match SNe 2017cbv and 2018oh, both of which had early light-curve excesses. We compare these models to our observed color curves using the same prodecure as above and find that the best-matched model requires a shell of 0.02 M_☉ of ⁵⁶Ni centered at 1.35 M_☉ with a thickness of 0.18 M_☉ in mass coordinates, offset by 0.5 days, which has ${\bar{\chi }}^{2}=116$ with N = 25. Again we see a reasonable match to the observed U − B and B − V color curves, but this time the bump is not strong enough to match the observed light curve (Figure 5, right).

For both the double-detonation and nickel-shell models, we were not able to find a better match by comparing the light curves directly, rather than the color curves. There are better matches to the light-curve peaks, but these do not correctly predict the early excess. Comparing these models to SN 2021aefx yielded the same problem (Hosseinzadeh et al. 2022).

Finally, we compared models of interaction with CSM from Piro & Morozova (2016) and models of pulsational delayed detonations from Dessart et al. (2014) to our observed light curves and color curves using the procedure above. However, the models in these sets are significantly poorer matches to our data than the three models discussed previously, so we do not discuss them further.

The features in our earliest spectra indicate very high velocities for the material ejected in the explosion. For example, the blue edge of the Si ii line at 636 nm in our first spectrum implies a maximum velocity of ∼42,000 km s⁻¹. The asymmetry of this feature in spectra more than 1 week before maximum light likely indicates two or more velocity components, as well as blending with the Si ii line at 597 nm. The broad, blended, hydrogen-free features of very young SNe Ia, in combination with the scarcity of equally early spectra to compare to in the literature, have led to mistaken Type Ic or uncertain Type I classifications in this event and others (Bostroem et al. 2021; Zhai et al. 2023).

We measured the expansion velocity and equivalent width of the Si ii line at 636 nm in each of our spectra by fitting a single Gaussian line profile in absorption on a linear continuum using an MCMC routine implemented in the emcee package (Foreman-Mackey et al. 2013). This is not always a perfect representation of the complete line profile, but it is sufficient to measure the velocity of the absorption minimum. In cases where C ii is also visible at 658 nm, we simultaneously fit a second Gaussian in absorption. The velocities are derived from the center of the Gaussian component of the model, and the equivalent widths come from integrating the data, after subtracting the model of the other line when necessary, divided by the model continuum. These quantities are also reported in Table 1.

Absorption due to C ii is clearly visible in our first five spectra (>14 days before peak; <3 days after explosion), with photospheric velocities around 21,000 km s⁻¹ and equivalent widths around 3.3 nm. Carbon is a signature of unburned material, which is not predicted by all explosion models, and has been previously observed in a significant fraction of early SN Ia spectra (e.g., Parrent et al. 2011; Thomas et al. 2011; Folatelli et al. 2012; Silverman & Filippenko 2012; Maguire et al. 2014).

Another remarkable feature of our measurements is that the equivalent width of Si ii initially increases from 9.8 ± 0.9 nm in our first spectrum to 14.1 ± 1.2 nm a few days later, followed by a decline until around maximum light. The shallow absorption lines in our earliest spectra suggest dilution of the lines by an additional continuum component of the emission. Marion et al. (2016) similarly noted shallow spectral features in spectra of SN 2012cg, which also showed an early light-curve excess, at these same phases. It is unclear whether this is normal behavior; for example, the earliest equivalent width measurements in the sample of Silverman & Filippenko (2012) are ∼13 days before maximum light. The extra continuum could come from several of the processes we discuss above: collision with a binary companion, extra heating from surface nickel, or interaction with CSM.

Ashall et al. (2022) postulated that the early U/u-band excess in SN 2021aefx was caused purely by the strong Ca ii H and K absorption feature traversing the U/u transmission function as the photosphere decelerates with time. However, the earliest observed velocity in SN 2023bee is ∼25,000 km s⁻¹, which is near the minimum requirement for a U-band bump according to Ashall et al. (2022). This mechanism is not able to explain excesses in most other bands, including the clear deviation of the unfiltered light curve from a power law, as well as the weaker excesses in all optical bands implied by our fitting.

To demonstrate the effect of SN 2023bee's lower early velocity on the early U − B color in the "speed bump" paradigm of Ashall et al. (2022), we blueshift our SALT spectrum taken 15.2 days before peak (the first-night spectrum that extends farthest into the blue) from a photopheric velocity of 24,000 to 30,000 km s⁻¹ to simulate the earliest velocity of SN 2021aefx. We then measure synthetic photometry in U and B on the original and blueshifted spectra. In U, where our spectrum does not cover the full passband, we only integrate over the wavelengths covered by the spectrum. We find that the higher-velocity spectrum is bluer by Δ(U − B) = − 0.2 mag, which is approximately how much bluer SN 2021aefx is than SN 2023bee. Considered alongside the evidence presented above, we conclude that, while this velocity effect may cause SNe Ia to be bluer during the first few days after explosion, particularly in the U and u bands, it cannot be the sole explanation for the early light-curve behavior of either of these events.

We can independently assess possible single-degenerate progenitors by exploring what properties of wind-like CSM are allowed by our radio upper limits (Section 2.5). Following the methodology of Chomiuk et al. (2012, 2016), we assume SN 2023bee was surrounded by the Chevalier (1982) model of CSM, produced by steady mass loss from the progenitor, i.e., ${\rho }_{\mathrm{csm}}=\dot{M}/4\pi {r}^{2}{v}_{{\rm{w}}}$, where ρ_csm is the CSM density, $\dot{M}$ is the mass-loss rate from the progenitor, r is the distance from the progenitor, and v_w is the wind velocity. We fix the energy fraction in relativistic electrons to _e = 0.1. We note that Chomiuk et al. (2012) assumes synchrotron self-absorption at low frequencies, but at early times such as our VLA observation, free–free absorption can also become important for high-density winds. We modify the optically thin luminosity from Chomiuk et al. (2012) by a factor ${{\rm{e}}}^{-{\tau }_{\mathrm{ff}}}$ (where τ_ff is the free–free optical depth) to account for free–free absorption. We use the definition of τ_ff given in Fransson et al. (1996) as

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=17.8\ {\nu }_{6}^{-2}{\dot{M}}_{-5}^{2}{v}_{{\rm{w}},1}^{-2}{T}_{5}^{-3/2}{({v}_{4}{t}_{10})}^{-3},\end{eqnarray} \tag{ 1 }$$

where ${\dot{M}}_{-5}=\dot{M}/{10}^{-5}\ {M}_{\odot }{\mathrm{yr}}^{-1}$ is the mass-loss rate, v_w,1 = v_w/10 km s⁻¹ is the wind speed, T₅ = T/10⁵ K is the wind temperature, v₄ = v/10⁴ km s⁻¹ is the ejecta velocity, t₁₀ = t/10 days is the SN age, and ν₆ = ν/6 GHz is our observing frequency. We assume T₅ = 1 for our calculations. ²⁴ The resulting light curves for different mass-loss rates are shown in Figure 3.

Assuming a standard SN Ia explosion with 10⁵¹ erg kinetic energy and 1.4 M_☉ ejecta mass, we obtain a mass-loss rate upper limit of $\dot{M}/{v}_{{\rm{w}}}\lt (6.4-30.5)\,\times {10}^{-10}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\ {\mathrm{km}}^{-1}\ {\rm{s}}$. The range of mass-loss rates reflects the uncertainty in the parameter _B , the fraction of shock energy shared by the amplified magnetic field, with typical values in the range 0.01–0.1 for SNe (Chomiuk et al. 2012). These limits are also consistent with the reported upper limits from ATCA (Leung et al. 2023), particularly the 5.5 GHz (shown in Figure 3) and 9 GHz limits, ²⁵ which, though shallower than our 6 GHz VLA data, are similarly constraining, because they were taken a few days earlier when the light curve would have been peaking.

These limits are compared with the mass-loss-rate parameter space of single-degenerate models, as defined by Chomiuk et al. (2012), in Figure 3. We find that our limits are deep enough to rule out the majority of red-giant companions (symbiotic systems), characterized by slow winds of 10–100 km s⁻¹ and mass-loss rates of 10⁻⁶–10⁻⁸ M_⊙ yr⁻¹ (Seaquist & Taylor 1990). Symbiotic systems have also been ruled out for the majority of SNe Ia based on radio upper limits (Chomiuk et al. 2012; Horesh et al. 2012; Pérez-Torres et al. 2014; Chomiuk et al. 2016; Lundqvist et al. 2020; Pellegrino et al. 2020; Burke et al. 2021). Many models involving main-sequence companions for SN 2023bee, however, are still allowed within our limits (see the colored regions in Figure 3). Weak constraints are obtained for the channel of fast, optically thick winds from white dwarfs accreting at high rates, while more quiescent mass-loss and nova-like densities are still allowed by our observations.

We briefly note that, for mass-loss rates >3 × 10⁻³ M_☉ yr⁻¹ (assuming v_w = 100 km s⁻¹), free–free absorption of radio emission will still yield nondetection at our 85 day ATCA observations at 5.5 and 9 GHz (for the parameters in Equation (1)). However, such high mass-loss rates in thermonuclear SNe have only been observed in a rare variety of Type Ia CSM, which shows signatures of dense CSM interaction, such as broadened emission lines of hydrogen or helium in the spectra (Silverman et al. 2013). No such features are seen in the spectra of SN 2023bee.

Our radio upper limits therefore make it unlikely that a red-giant companion was part of SN 2023bee's progenitor system, unless the ratios of energies in the magnetic fields and relativistic electrons to the shock energy are much smaller than typically assumed (i.e., _B , _e ≪ 1%–10%). A red-giant companion is also effectively ruled out by our light-curve fitting in Section 3. However, the main-sequence companion preferred by our light-curve fitting is still permitted by our radio upper limits.

Taken together, our observations above suggest that the emission from SN 2023bee during the first few days after explosion consists of typical SN Ia spectral features plus an additional hot continuum component. This causes an excess in all bands but preferentially at blue to UV wavelengths. Of the models we consider, the companion-shocking model of Kasen (2010) best reproduces the strength, duration, and color of the excess in the UVW1 through the i bands. Other models, including those involving a double detonation or a shell of ⁵⁶Ni near the surface of the ejecta, are able to produce an excess of approximately the same duration but do not match in either strength or color. The worse match of these latter models might be attributed to their relatively small grid of parameters, while the analytic companion-shocking model is more flexible. Larger model grids or analytic approximations could allow for fairer comparisons in the future.

One weakness of the Kasen (2010) model, already noted in previous work (Hosseinzadeh et al. 2017, 2022), is the assumption that the shock component of the emission does not suffer from the same line blanketing in the UV as the emission from the rest of the ejecta. This leads to an overprediction in the farther-UV bands (UVM2 and UVW2). The severity of this overpredicton likely depends on the details of the spectrum, including line strengths and velocities, which in turn depend on the conditions in the ejecta and possibly the viewing angle. Future modeling work should focus on a more realistic spectral energy distribution for the shock component with the goal of a fairer comparison with the data.

In this single-degenerate scenario, the ejecta are expected to strip material off the companion star during the collision, and this material could emit at late phases when the ejecta are mostly transparent (Botyánszki et al. 2018; Dessart et al. 2020). SN 2023bee has not yet reached this phase, but the search for nebular emission will be an important test of the single-degenerate model. Previous searches for hydrogen and helium in nebular spectra of SNe Ia have mostly yielded null results (Leonard 2007; Graham et al. 2017; Sand et al. 2019; Tucker et al. 2020), even in cases where there was an early light-curve excess (Sand et al. 2018; Dimitriadis et al. 2019b; Tucker et al. 2019; Siebert et al. 2020; Tucker et al. 2021; Hosseinzadeh et al. 2022). However, translating line luminosity limits to limits on the stripped mass depends on complex and poorly understood physical processes.

This Letter demonstrates the power of using very high-cadence, multiband photometry of young, nearby SNe Ia to constrain their progenitor systems, which is only possible with specially designed robotic facilities like the DLT40 Survey, Las Cumbres Observatory, and Swift. As the number of well-observed events grows, increasingly high-quality data sets are beginning to reveal the weaknesses in existing models of SN Ia progenitors and explosion physics. Close cooperation between observers and modelers is required to narrow this gap.

This work makes use of data from the Las Cumbres Observatory telescope network. The SALT data presented here were taken via the Rutgers University program 2022-1-MLT-004 (PI: S.W.J.) and the South African program 2021-2-LSP-001 (PI: D.A.H.B.). Polish participation in SALT is funded by grant No. MEiN nr2021/WK/01. The Australia Telescope Compact Array is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Gomeroi people as the traditional owners of the Observatory site. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Time domain research by G.H, D.J.S, and the University of Arizona team is supported by NSF grants AST-1821987, 1813466, 1908972, and 2108032, and by the Heising-Simons Foundation under grant #2020–1864. J.E.A. is supported by the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the United States of America. This publication was made possible through the support of an LSSTC Catalyst Fellowship to K.A.B., funded through Grant 62192 from the John Templeton Foundation to LSST Corporation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of LSSTC or the John Templeton Foundation. Research by S.V., Y.D., E.H., and N.M. is supported by NSF grant AST-2008108. The Las Cumbres Observatory group is supported by NSF grants AST-1911225 and 1911151. L.A.K. acknowledges support by NASA FINESST fellowship 80NSSC22K1599. D.A.H.B. acknowledges research support from the National Research Foundation of South Africa. L.C. is grateful for support from NSF grant AST-2107070.

Facilities: ADS - , ATCA - Australia Telescope Compact Array, Bok (B&C) - Steward Observatory 2.3 meter Bok Telescope, CTIO:PROMPT - Panchromatic Robotic Optical Monitoring and Polarimetry Telescopes at Cerro Tololo Inter-American Observatory, LCOGT (Sinistro - , FLOYDS) - , Meckering:PROMPT - , MMT (Binospec) - MMT at Fred Lawrence Whipple Observatory, NED - , SALT (RSS) - Southern African Large Telescope, SOAR (Goodman) - The Southern Astrophysical Research Telescope, Swift (UVOT) - Swift Gamma-Ray Burst Mission, VLA - Very Large Array

Software: Astrometry.net (Lang et al. 2010), Astropy (Astropy Collaboration et al. 2022), BANZAI (McCully et al. 2018), CASA (McMullin et al. 2007), corner (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013), extinction (Barbary 2016), FLOYDS pipeline (Valenti et al. 2014), Goodman pipeline (Torres et al. 2017), HEAsoft (NASA HEASARC 2014), IRAF (National Optical Astronomy Observatories 1999), lcogtsnpipe (Valenti et al. 2016), Light Curve Fitting (Hosseinzadeh et al. 2023a), Matplotlib (Hunter 2007), NumPy (Oliphant 2006), Photutils (Bradley et al. 2022), PyRAF (Science Software Branch at STScI 2012), PySALT (Crawford et al. 2010), reproject (Robitaille et al. 2023), SNCosmo (Barbary et al. 2022).