Constraints on the Progenitor of SN 2023bee from its Nebular-Phase Spectra

Corresponding author: Sophie von Coelln

[email protected] ferred velocities validated against NIR spectrum at the similar phase. The resulting low Ni-to-Fe mass ratio, $m_{\mathrm{Ni}}/m_{\mathrm{Fe}} \lesssim0.06$, is broadly consistent with a sub-$M_{\mathrm{Ch}}$ double-detonation progenitor system.

Introduction

Type Ia supernovae (SNe Ia) are generally believed to originate from the thermonuclear explosion of a white dwarf (WD) star, yet the nature of their progenitor systems remains uncertain. In particular, the properties of the progenitor system specifically the type of companion star, as well as the physical conditions leading to explosion, are not yet fully established. This uncertainty is particularly important given the use of SNe Ia as standardized candles for measuring the expansion history of the Universe [99, 87]. Different progenitor and explosion scenarios are expected to produce different early-time signatures, including excess emission shortly after explosion, but uncertainties in theoretical modeling and limitations in early-time follow-ups often make it hard to exclusively connect them with certain progenitor systems. Nebular-phase spectroscopy provides complementary constraints on the inner ejecta and can help distinguish among otherwise degenerate scenarios.

The complex physics behind SN Ia explosions has been investigated for decades [e.g. 52, 71, 101, 70, 102]. The two main channels for SN Ia progenitor systems are single degenerate (SD) and double degenerate (DD) systems. In the SD scenario, the WD explodes due to a thermonuclear runaway triggered as the WD accretes mass from a nondegenerate companion star near the Chandrasekhar mass ($M_{\mathrm{Ch}}$; [126]). In the DD scenario, mass transfer or a merger event between two WDs results in a thermonuclear runaway, leading to a standard or peculiar supernova explosion [55, 31, 45, 102]. These broad progenitor classes encompass a range of binary configurations and evolutionary pathways that may produce both normal and peculiar SNe Ia.

A separate but related question concerns the mass and ignition conditions of the exploding WD. Near-Chandrasekhar-mass ($M_{\mathrm{Ch}}$) explosions occur when the WD grows to approximately $M_{\mathrm{Ch}}$ before thermonuclear runaway [e.g. 109, 19], whereas sub-$M_{\mathrm{Ch}}$ explosions occur at substantially lower masses. Sub-$M_{\mathrm{Ch}}$ explosions could be triggered by double detonation [e.g. 12, 111, 91, 86, 15] or a merger [83, 84, 103]. Other proposed channels include the core-degenerate (CD) scenario [54, 105, 59, 123], distributed burning realized in turbulently driven detonation regime [127, 93, 18, 129], and all versions of collisions such as head-on or impact parameter between double WDs [8, 96, 63, 30, 85, 41]. More generally, the thermonuclear burning may proceed through deflagration, detonation, or transitions between burning regimes, further broadening the range of possible SN Ia explosion outcomes.

One potential way to distinguish between SD and DD scenarios is the early light curve: we expect the early light curve of SN Ia to rise roughly following a power law, with the power-law index close to 2 in the absence of additional emission [2, 47]. Several mechanisms can produce an early excess, or “bump” in extreme cases, including ejecta–companion interaction in the SD scenario [58]; shallow or asymmetric ${}^{56}\mathrm{Ni}$ distributions [89, 73, 72]; circumstellar interaction [88, 80]; and radioactive products of He-shell detonations in double-detonation models [91]. Distinguishing among these possibilities requires densely sampled observations within the first hours to days after explosion. Such data are difficult to obtain because most SNe are discovered after the earliest phases, while rapid follow-up is not always available. In addition, uncertainties in explosion and radiative-transfer models, together with degeneracies among model parameters and viewing angle, limit the extent to which even well-observed early excesses can be uniquely interpreted.

In addition to the observational features like the early light curve, nebular phase spectroscopy provides another window into the progenitor nature of SNe Ia. In an SD scenario, where the SN ejecta shock the nondegenerate companion, they will also strip material from the surface of the companion, which will be evident in spectra taken during the nebular phase of the SN [78, 67, 14, 16]. The stripped material will produce strong and narrow emission lines on top of an otherwise standard SN Ia spectrum. In particular, we expect to see some or all of H$\alpha$, H$\beta$, He I $\lambda5875$, and He I $\lambda6678$, varying based on the type of companion star and the characteristics of the binary system.

Moreover, the density attained during thermonuclear burning affects the degree of electron capture and therefore the production of stable neutron-rich iron-group isotopes, providing a potential diagnostic of the SN Ia explosion mechanism [46, 77]. In particular, nebular [Fe, II] and [Ni, II] emission can be used to constrain the stable Ni-to-Fe abundance ratio, which is expected to differ between near-$M_{\mathrm{Ch}}$ and sub-$M_{\mathrm{Ch}}$ explosion scenarios [34, 113]. We discuss these diagnostics and their associated uncertainties in detail in § 4.

A number of SNe Ia have been found with early flux excess, including SNe 2017cbv, 2018oh, 2021aefx [49, 28, 110, 4] (Li et al. 2019). These SNe range from normal to overluminous, and their flux excess is normally blue and bright in UV with a few red exceptions (such as SN 2018aoz; Ni et al. 2023). The shape of 2018oh’s early light curve can be fit by the SD companion shocking model, the shallow ${}^{56}\mathrm{Ni}$ model, or the sub-M$_{\mathrm{Ch}}$ double detonation model, but none of these models can explain both the blue early excess and the lack of narrow emission lines in nebular spectroscopy. SNe 2017cbv and 2021aefx have early UV luminosity measurements from Swift, which are too low to align with the bright UV bump expected from the SD companion-shocking model. In addition, SD models have been largely excluded for SNe 2017cbv and 2018oh by multiple studies, due to insufficient hydrogen and helium lines in their nebular spectra ([104, 29, 118]). Alternative explanations also face observational challenges, models containing an outer ${}^{56}\mathrm{Ni}$ shell can reproduce the early optical excess but fail to match the color and spectral evolution, predicting excessive line blanketing and an unobserved absorption feature around 3700 – 4000 Å ([72]). Interaction with H/He-poor CSM remains a possible explanation for SN 2017cbv, and interaction with disk-originated matter has been shown to reproduce the early rise of SN 2018oh (e.g., Hosseinzadeh et al. [49]). None of the current models are sufficient to describe potential progenitor systems for all SNe Ia with early excess flux.

In this paper, we examine the late-time spectra of SN 2023bee, a normal SN Ia observed by the Transiting Exoplanet Survey Satellite (TESS) discovered with a weak early excess. It was discovered 2 days post-explosion and also observed photometrically for 12 days pre-explosion with a 10-minute cadence, identified with . This extensive observation allows us to closely investigate the early flux and nebular spectrum of SN 2023bee and their effects on the progenitor system model for the SN. Hosseinzadeh et al. [51] presents an early spectral time series and radio observations of the SN, concluding that the companion-interaction model in Kasen [58] best aligns with the observed early excess flux, and noting that, under such a model, the ejecta should strip material from the companion, which could be observed in nebular-phase spectroscopy. Wang et al. [124] presents a more detailed examination of the photometric data, including the TESS photometry, finding the limitations of existing models in explaining the early excess of SN 2023bee across all UV-optical bands. This motivates additional late-time observations to better constrain the progenitor system of SN 2023bee. In § 2 we present the acquisition, reduction and calibration of our photometric and spectroscopic data. In § 3 we analyze the H/He features in the optical and NIR spectra, and place upper limits on the H/He-rich material stripped from a companion. In § 4 we estimate the iron-group abundance ratio to constrain the central density and the progenitor system. We present the implications for the progenitor sys-

tem in § 5. All the data presented in this paper will be made public via WIS eREP ¹ (Yaron & Gal-Yam [128]).

Data Acquisition, Reduction and Calibration

SN 2023bee was discovered at coordinates $\alpha= 08{:}56{:}11.63$ and $\delta= -03^\circ19' 32.06''$ (J2000.0). It peaked at 59992.58 MJD with an apparent magnitude of $m_B = 13.041 \pm0.001$ mag, and an absolute magnitude of $M_B = -19.6 \pm0.2$ mag (Wang et al. [124]). It is located in galaxy NGC 2708, an intermediate spiral galaxy at $z = 0.0067 \pm0.0005$, and is 93.97″ from the center of the galaxy. It has a distance modulus $\mu= 32.5 \pm0.2$ mag (Pisano et al. [90]). We expect negligible host extinction given SN 2023bee’s large projected offset from host nucleus. In this paper, we use $E(B-V)_{\mathrm{MW}} = 0.145$ for Milky Way extinction based on the map given in Schlafly & Finkbeiner [107], and dust extinction law from Fitzpatrick [32].

Late-time Photometry

We obtain the late-time photometry of SN 2023bee with the Nickel 1.0-m telescope at the Lick Observatory in $BVri$ bands, the Swope 1.0-m telescope at Las Campanas Observatory in $uBVgri$ bands, and the Pan-STARRS1 1.8 m telescope at Haleakala on Maui, Hawai’i in $w$ bands (Chambers et al. [20]), as part of the Young Supernova Experiment (YSE; Jones et al. [56]). We also include the publicly available data from the Zwicky Transient Facility (ZTF; Bellm et al. [7]; Dekany et al. [23]) in $gr$ bands. Standard reductions for the DECam and PS1 images are performed by the NOIRLab community pipeline (Valdes et al. [121]) and the PS1 Image Processing Pipeline (IPP), respectively (Magnier et al. 2020a,b,c; Waters et al. 2020). These images are then taken as input to the photpipe pipeline (Rest et al. [97]; Rest et al. [98]), which re-determines the zero points by comparing DoPHOT PSF photometry from each image to the Pan-STARRS DR1 catalog (Flewelling et al. [33]), convolves and subtracts a template image from the survey image, and performs forced photometry on the resulting difference images to create transient light curves. As described in Kilpatrick et al. [60], all Swope image processing and optical photometry was performed using photpipe (Rest et al. [97]), flat-fielding, image stitching, and photometric calibration. $B\!Vri$ photometry were calibrated using standard sources from the Pan-STARRS DR1 catalog (Flewelling et al. [33]), while the $u$-band data were calibrated using SkyMapper $u$-band standards (Onken et al. [82]), transformed into the Swope natural system (Krisciunas et al. [62]) with the Supercal method (Scolnic et al. [108]). Nickel $BVri$ imaging was reduced with photpipe (Rest et al. [97]; Kilpatrick et al. [60]) following the procedures described for Nickel and LCO data in Gagliano et al. [40] and Kilpatrick et al. [61]. We applied bias and flat-field frames obtained on the same night and in the same instrumental configuration, aligned each frame to 2MASS astrometric standards [116], and regridded the images with SWarp [9] to remove geometric distortion. PSF photometry was performed with a custom version of DoPHOT [106]. The $gri$ photometry was calibrated to Pan-STARRS DR1 standards [33] in the same field, and Johnson $BV$ magnitudes were obtained by transforming those PS1 standards following Tonry et al. [117].

Optical and NIR Spectroscopy

We report two collected late-time optical spectra of SN 2023bee in this paper. One optical spectrum was taken by the Kast spectrograph on the Lick 3 m (Shane) telescope on 2023 May 11 (MJD 60075), $\sim+82$ days after the peak in rest frame. Another optical spectrum was taken by Keck telescope with LRIS spectrograph on Feb 14 2024 (MJD 60354), $\sim+359$ days after the peak in nebular phase. To reduce the Kast and LRIS data, we used the UCSC Spectral Pipeline² [112], a custom data-reduction pipeline based on procedures outlined by Foley et al. [35], Silverman et al. [114], and references therein. The two-dimensional (2D) spectra were bias-corrected, flat-field corrected, adjusted for varying gains across different chips and amplifiers, and trimmed. One-dimensional spectra were extracted using the optimal algorithm [48]. The spectra were wavelength-calibrated using internal comparison-lamp spectra with linear shifts applied by cross-correlating the observed night-sky lines in each spectrum to a master night-sky spectrum. Flux calibration and telluric correction were performed using standard stars at a similar airmass to that of the science exposures. We combine the red and blue sides by scaling one spectrum to match the flux of the other in the overlap region and use their error spectra to correctly weight the spectra when combining. More details of this process are discussed elsewhere [35, 114, 112]. Two NIR spectra were taken by Magellan telescope with the folded port infrared echelle (FIRE; [115]). They were reduced and extracted with PyPeIt (Prochaska et al. 2020; Prochaska et al. 2020). We chose to rebin the LRIS and FIRE spectrum to 3.5Å and 1.67Å respectively to approximately match their actual spectroscopic resolution power. The Kast spectrum was not rebinned as it only serve as photometric calibrator in later analysis.

The spectrum at +359 day for SN 2023bee is shown in light gray in Figure 2; overplotted in black is the smoothed spectrum, and absolute flux calibration is discussed in § 2.3. Also shown in blue are the normalized, redshift-corrected, and similarly smoothed nebular phase spectra of three SNe Ia with early excess, i.e. SN 2017cbv, SN 2018oh and SN 2021aefx, at +286, +264 and +303 days respectively. The LRIS optical spectrum taken on +359 days shows blends of broad features from ionized Fe, Co and Ni, includ-

MJD (day)	Phase$^{a}$ (day)	Filter	$m$ (mag)	$\sigma_m$ (mag)	Telescope/Instrument
59976.998	-15.478	o	16.793	0.030	ATLAS/HKO
59977.018	-15.458	o	16.764	0.017	ATLAS/HKO
59978.552	-13.935	g	15.510	0.010	Pan-STARRS1/GPC1
59978.552	-13.935	r	15.649	0.011	Pan-STARRS1/GPC1
59978.553	-13.934	i	16.017	0.014	Pan-STARRS1/GPC1
59978.554	-13.933	z	15.848	0.014	Pan-STARRS1/GPC1
59978.555	-13.932	y	16.105	0.034	Pan-STARRS1/GPC1
59979.294	-13.198	g	15.242	0.024	P48/ZTF-Cam
59979.401	-13.091	r	15.309	0.025	P48/ZTF-Cam
59981.430	-11.076	g	14.250	0.005	Pan-STARRS1/GPC1
59981.431	-11.075	r	14.406	0.005	Pan-STARRS1/GPC1
59981.431	-11.075	i	14.829	0.005	Pan-STARRS1/GPC1
59981.432	-11.074	z	14.812	0.007	Pan-STARRS1/GPC1
59981.433	-11.073	y	15.055	0.017	Pan-STARRS1/GPC1
59983.204	-9.314	g	13.868	0.016	P48/ZTF-Cam
59983.337	-9.181	r	14.000	0.015	P48/ZTF-Cam
59983.397	-9.122	g	13.704	0.002	Pan-STARRS1/GPC1
59983.398	-9.121	r	13.907	0.004	Pan-STARRS1/GPC1
59983.398	-9.121	i	14.254	0.005	Pan-STARRS1/GPC1
59983.399	-9.120	z	14.360	0.008	Pan-STARRS1/GPC1

ing the strong [Fe III] emission around 4700 Å and [Co III] around 5900 Å. Another noticeable feature is the blended region around 7200 Å composed of [Fe II] and [Ni II] multiplets, and potentially existing [Co III] and [Ca II] features. As shown in Figure 2, SN 2023bee show no detectable narrow H/He features that are expected in the SD scenario, consistent with other well-studied SNe Ia exhibiting early excesses, as illustrated in the bottom residual panel.

Similar to optical spectra, the +101 day NIR spectrum of SN 2023bee shows strong broad iron-group line features including [Co II] 10190 Å, [Fe II] 12570 Å and [Co II] 15470 Å. At +333 day, the iron-group line features become more prominent relative to weaker continuum, with broad [Fe II] 16440 Å feature also appearing. On the other hand, no significant H/He feature can be seen on these nebular phase spectra of SN 2023bee. The analysis on the line features around 7200 Å and in NIR will be discussed in § 4.

Spectroscopic flux calibration

To obtain an accurate line intensity measurement for further physical analysis, the spectra need to be flux calibrated.

The spectroscopic synthetic magnitude in each filter was determined with the PYPHOT package from the fluxes and wavelengths of the two optical spectra. The photometric magnitude was then found via linearly interpolating late-time photometry from the Swope $gri$ filters to infer the accurate magnitudes of SN 2023bee in each band on the date the spectrum was taken, as shown in the Figure 1. We re-port the interpolated magnitudes to be $m_g = 20.21$, 22.14 and 21.23 mag in the $gri$ bands respectively at +359 day, in good agreement with neighboring Swope photometric points with $m_g = 20.208$ mag and $m_i = 21.2$ mag at +358 day, $m_r = 22.089$ mag at +362 day. Computing the difference between the photometric and spectroscopic synthetic magnitude in each band yields -37.306, -37.231, and -37.188 mags in $g$, $r$, and $i$ bands, respectively. The $\sim0.1$ mag difference between $g$ and $r/i$-band offsets may arise from a relative flux-calibration mismatch between the blue and red arms of LRIS, and we therefore include this difference as a component of the photometric calibration uncertainty. From there, we normalize the spectra with the average difference between the photometric and synthetic magnitude in $gri$-bands, and then correct for the distance modulus. Finally, the rest-frame spectrum was smoothed with a third-order Savitzky-Golay filter with window size of 100 Å. The small window size means this smoothed out noise while maintaining more significant features, allowing us to better identify any noticeable emission lines that may come from the remaining H and He from the progenitor. We estimate that the flux conversion scale of the +359 day optical spectrum has a $\sim4.7\%$ uncertainty associated with it given the uncertainties and scatter in photometric measurements.

Afterwards, we first flux-calibrate the optical spectrum on +80 day with $gri$ photometry in a similar way, and then flux-calibrate the adjacent +101 day NIR spectra based on the overlapping wavelength region of 8350–9200 Å. The 9200–10000 Å region is affected by telluric contamination, particularly in the Lick optical spectrum obtained on 2023 May 11 due to high-humidity conditions of that night, and was therefore excluded from the flux calibration. The magnitude offsets in the optical–NIR overlap region at different epochs were approximated using the corresponding changes in the $i$-band magnitudes. We caution, however, that between +80 and +101 day the line features such as Ca II triplets may still undergo non-negligible evolution ([39, 120]), which could introduce systematic biases in the absolute flux calibration. We adopt a linear flux-calibration relation and show the calibrated NIR spectra in physical units in Figure 3. On the other hand, owing to the large uncertainties in the overlapping region of the +333 day NIR spectrum below 10000 Å, we were unable to establish a reliable relative flux calibration with the +359 day optical spectrum. We therefore exclude this spectrum from subsequent analyses that rely on an accurate absolute flux scale.

Narrow H/He Features and non-degenerate companion

Obs Date (UT)	MJD	$a_{\mathrm{Phase}}$ (Days)	Telescope/Instrument
2023-05-11	60075.2	82.0	Shane/Kast
2023-06-01	60096.98	103.7	Magellan/FIRE
2024-01-19	60329.12	334.3	Magellan/FIRE
2024-02-14	60354.4	359.4	Keck/LRIS

Figure 4 shows zoom-ins for the +359 day LRIS spectrum for the H$\alpha$, H$\beta$, He I $\lambda5875$ and He I $\lambda6678$ lines. No prominent emission features are found at the expected wavelengths of relevant H and He features in optical and NIR in SN 2023bee. Therefore, we follow the methodology of Dimitriadis et al. [29] to place upper limits on H$\alpha$, H$\beta$, Paschen lines, and He I emission and thereby constrain the amount of stripped material. The LRIS spectrum is rebinned to 3.5Å to approximate LRIS spectral resolution. We estimate the spectral uncertainties from the residual obtained by subtracting a smoothed continuum of the observed spectrum and use them to derive quantitative upper limits on undetected line emission at these wavelengths. To determine an appropriate continuum while minimizing the influence of any potentially present weak, narrow emission lines, we apply a second-order Savitzky–Golay filter with a broad window of 260Å.

The LRIS optical spectrum obtained on +359 day were observed in red and blue filters separately, separated at around 5700Å. We therefore estimate the residual noise independently within the blue and red wavelength ranges. Specifically, we calculate the standard deviation of the residuals over 4160–5560Å for the blue spectrum and 5700–9200Å for the red spectrum. These intervals encompass the H$\alpha$, H$\beta$, Paschen lines, and He I $\lambda5875$ and He I $\lambda6678$ regions, and the uncertainty levels are relatively uniform within each of the blue or red spectral regions. We therefore adopt the global uncertainties measured in red and blue arms rather than the local uncertainties around the individual line features. To avoid biasing the noise estimate from any weak emission at the expected line positions, we exclude ±100Å around each transition. We obtain a $3\sigma$ flux-density upper limit of $8.3 \times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ for H$\beta$ in the blue filter. For the red filter, we find a corresponding $3\sigma$ upper limit of $2.0 \times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ for Hα, He I λ5875, and He I λ6678.

To visualize the limit, we add narrow Gaussian profiles on the optical spectrum at the rest wavelengths of Hα, Hβ, and He I λ5875 and He I λ6678 line regions. The widths of the Gaussian profiles are fixed at values corresponding to 1500 km s$^{-1}$ at the wavelengths of the respective emission lines, equivalent to $\sim33$ Å for Hα. The bottom panels of subplots in Figure 4 show the comparison between the continuum-subtracted residuals and those artificial narrow Gaussian profiles representative of the corresponding flux limits. The comparison confirms that no significant narrow emission features are present at the wavelengths of interest. The resulting flux upper limits are 5.4, 1.6, 1.8, and $1.8 \times10^{37}$ erg s$^{-1}$ for Hβ, He I λ5875, Hα, and He I λ6678. We caution, however, that the Hβ limit is conservative. As shown in the upper right panel of Figure 4, the broad Fe features in the vicinity of the Hβ region introduce additional residuals relative to the smoothed continuum. Consequently, the estimated uncertainty near Hβ line is an overestimate of the intrinsic spectral noise.

Botyánszki et al. [17] derived the empirical relation between the line intensity and amount of stripped mass with an analytical fit to the multidimensional non-LTE radiative-transfer calculations basd on ejecta–companion interaction model in the Roche-lobe-filling scenario (Boehner et al. [14]). With (1) in Botyánszki et al. [17],

where $M_{1} = \log_{10}(M_{\mathrm{st}}/M_{\odot})$ and $[L_{\mathrm{H}\alpha}]$ is in unit of erg s$^{-1}$, we estimate the upper limit of stripped H to be $M_{\mathrm{H}}^{\mathrm{stripped}} \lesssim5.0 \times10^{-4}M_{\odot}$. Botyánszki et al. [17] also estimated that the same equation will overestimate the mass of stripped He by a factor of four in a He-only model. Therefore, we place an upper limit of the stripped He as $M_{\mathrm{He}}^{\mathrm{stripped}} \lesssim1.2 \times10^{-4}M_{\odot}$.

Dessart et al. [26] pointed out that, if the stripped material is mixed within the ejecta when the spectra were taken, the optical depth will largely suppress H I features, making Hβ almost vanish before +300 days and decreasing Hα intensity by a factor of 3 – 5. Using the polynomial index fitted for three different simulation models at +300 days (see Appendix A of Dessart et al. [26]), we estimated an upper limit $M_{\mathrm{H}}^{\mathrm{stripped}} \lesssim5.0 \times10^{-4}M_{\odot}$ with the Hα feature in these scenarios. A similar effect also influences the He I features in optical, although, to our knowledge, no numerical models have yet quantified its magnitude in the case of He-rich ejecta.

Additionally, Dessart et al. [26] simulated nebular phase NIR spectra of SNe Ia in the H-rich companion scenario, including emission from stripped H/He-rich material, and predicted observable strong He I 10830Å and Paγ. Figure 5 compare the FIRE NIR spectrum at +101 day with the corresponding +100 day models that most closely match the observational detection limits in the He I 10830 Å and Pa$\gamma$ line regions. Following a similar approach, we derive an upper limit of $2.3 \times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ for Pa$\gamma$ and He I $\lambda$10830. In this comparison, the contribution from stripped materials as the excess flux relative to a baseline model with a minimal amount of stripped H-rich material, $M_{\mathrm{st}} = 1 \times10^{-5}M_{\odot}$. The +101 day NIR spectrum yields an upper limit of $M_{\mathrm{st}} \lesssim0.01M_{\odot}$ for stripped H-rich material, which is less constraining than the limit derived from the +359 day optical spectrum. The constraint from the FIRE NIR spectrum at +333 days is even weaker given its low S/N and large uncertainties in absolute flux calibration.

Iron-group Line Diagnostics

As discussed in § 1, the relative abundances of iron-group isotopes in the innermost ejecta are sensitive to the neutron excess and density attained during thermonuclear burning, and therefore provide important diagnostics of the explosion mechanism ([46, 10]). Near-$M_{\mathrm{Ch}}$ models with steady accretion can reach central densities of $\sim10^{9}$ g cm$^{-3}$, substantially higher than the $\sim10^{7}$ g cm$^{-3}$ typical of sub-$M_{\mathrm{Ch}}$ models ([100, 83, 1]). The higher densities enhance electron capture and the production of stable neutron-rich iron-group isotopes, such as ${}^{54}\mathrm{Fe}$ and ${}^{58}\mathrm{Ni}$, relative to radioactive ${}^{56}\mathrm{Ni}$ and its decay products ([46, 77]). The [Ni II]/[Fe II] flux ratio has therefore been proposed as a probe of neutronization and, in turn, of the distinction between near-$M_{\mathrm{Ch}}$ and sub-$M_{\mathrm{Ch}}$ explosions ([34]).

This diagnostic becomes accessible at nebular phases, $\gtrsim200$ days after maximum light, when the ejecta are sufficiently optically thin to reveal emission from the iron-group-rich inner layers. By these epochs, the short-lived radioactive isotopes ${}^{56,57}\mathrm{Ni}$ have decayed to ${}^{56,57}\mathrm{Co}$ and to ${}^{56,57}\mathrm{Fe}$ as:

At $\gtrsim300$ days, [Ni II] emission is expected to arise predominantly from stable neutron-rich Ni isotopes, particularly ${}^{58}\mathrm{Ni}$, with a smaller contribution from ${}^{60}\mathrm{Ni}$. In contrast, [Fe II] emission is dominated by the decay products of radioactive ${}^{56}\mathrm{Ni}$, with additional contributions from stable Fe isotopes and other radioactive decay chains ([42, 43, 77, 34, 112, 109]). Because the production of stable neutron-rich iron-group isotopes depends strongly on the density and electron fraction of the burning material, near-$M_{\mathrm{Ch}}$ and sub-$M_{\mathrm{Ch}}$ models generally predict different Ni/Fe abundance ratios. The [Fe II], [Ni II], and [Ca II] complex around 7300 Å can therefore constrain the stable Ni-to-Fe abundance ratio, with comparatively high Ni/Fe ratios generally associated with the enhanced neutronization expected in near-$M_{\mathrm{Ch}}$ explosions, whereas lower ratios are more readily reproduced by sub-$M_{\mathrm{Ch}}$ models ([34, 113]). Complementary NIR transitions, particularly [Ni II] $\lambda1.939\,\mu\mathrm{m}$ and [Fe II] lines such as $\lambda1.257\,\mu\mathrm{m}$, provide additional constraints on the stable Ni abundance and the physical conditions of the iron-rich ejecta ([34, 11, 10, 13]). In this work, we use the $\sim7200$ Å complex together with NIR iron-group lines as an independent diagnostic of the explosion scenario. The inferred Ni/Fe ratio must be interpreted with caution though as the line intensities also depend on the ionization state and the mixing of different isotopes in the ejecta ([10, 13]).

The line complex between 7000 -- 7600 Å is dominated by [Fe II] (7155, 7172, 7388, 7453 Å) and [Ni II] (7378, 7412 Å). [Co II] (6853, 6961, 7152 Å) also contribute to this region, particularly on its blue side, thereby affecting the inferred [Ni II] and [Fe II] line strengths at a non-negligible level ([34, 69]). Previous studies suggest that [Ca II] (7291, 7324 Å) is another component that could contribute to this region ([79, 38] (Błotyański & Kasen 2017)), though a number of empirical studies have put strong upper limits on its existence in other SNe Ia ([77, 42]). In the following analysis, we assume that the observed emission is dominated by these transitions and lines arising from the same ion share a common velocity offset and line width. The rest wavelength and Einstein coefficient of those lines are taken from the NIST Atomic Spectra Database.⁴

Figure 6 shows a zoom-in plot for the optical LRIS spectrum taken on +359 day in the $\sim7200$ Å complex region, normalized to the peak flux in this region. We fit this feature using several multi-component emission-line models: Fe+Ni, Fe+Co+Ni, and Fe+Co+Ni+Ca. For configurations including [Co III], we extend the blue boundary of the fitting region to $\sim6800$ Å in order to include the [Co III] features at $\lambda6854$ and $\lambda6961$, while the fitting range of Fe+Ni model is restricted 7000--7600 Å to exclude [Co III]. Our fitting procedure broadly follows the approach of Siebert et al. [113], with one important modification in the treatment of local continuum. In Siebert et al. [113], the continuum was estimated by linearly interpolating between spectral regions on either side of the emission complex. This approach can bias the inferred line fluxes when the adopted continuum regions themselves contain weak line emission, particularly on the blue side of the complex where [Co III] may contribute. We therefore adopt a constant continuum level across the fitting region and treat it as a free parameter.

!Best-fit for the 7300 Å region with [Fe II], [Ni II], and [Co III] bbox=0.08,0.08,0.48,0.37

The results for different fitting configurations are listed in Table 3. We compare the fitting results using both the reduced $\chi^2$ and the Bayesian information criterion (BIC). For Fe+Ni+Co and Fe+Ni+Co+Ca, we additionally report the reduced $\chi^{2}$ over both original fitting ranges and between 7000 − 7600 Å to enable a direct comparison with Fe+Ni model. The BIC is used primarily to compare Fe+Ni+Co and Fe+Ni+Co+Ca with different number of parameters, and the BIC for Fe+Ni is provided only for reference as the fitting range differs. In general, the Fe+Co+Ni model with a constant continuum shows the best performance while yielding consistent $v_{\mathrm{offset}}$ and line width among the iron-group species, and is therefore our preferred configuration. This result indicates that [Co III] makes a non-negligible contribution to the blended emission complex. When [Co III] is omitted, its contribution near the blue side of the complex is partially absorbed into the fitted [Fe II] component. This changes the inferred [Fe II] contribution at longer wavelengths and, in turn, reduces the flux attributed to [Ni II] λ7378. Consequently, models that neglect [Co, III] yield systematically lower inferred Ni/Fe ratios. Including [Ca II] produces somewhat more similar velocity offsets for the [Ni II] and [Fe II], but does not yield improvement in either $\chi^{2}$ or the BIC. Moreover, the inferred [Ca II] intensity remains below the $2\sigma$ level. We therefore find no statistically significant evidence for [Ca II] emission in this region, consistent with previous studies [77, 42]. We also test a linear continuum for each line-identification configuration. In all cases, the inferred continuum slopes are consistent with zero within the uncertainties. Allowing a linear continuum does not significantly improve either $\chi^{2}$ or the BIC, nor does it produce statistically significant changes in the best-fit [Fe II] or [Ni II] parameters.

Model	fitting range [Å]	Lines	Relative Amplitude	$v_{\mathrm{offset}}$ [km s$^{-1}$]	Width [km s$^{-1}$]	Continuum	$n_{\mathrm{NiII}}/n_{\mathrm{FeII}}$	$\chi^2_\nu{}^{a}$	BIC$^{b}$
Fe+Ni	[7000, 7600]	[Fe II]	$0.61 \pm0.01$	$875 \pm55$	$3490 \pm100$	$0.15 \pm0.01$	$0.0125 \pm0.009$	176	-500
Fe+Ni		[Ni II]	$0.063 \pm0.006$	$-50 \pm470$	$1880 \pm420$
Fe+Ni+Co	[6800, 7600]	[Fe II]	$0.60 \pm0.01$	$800 \pm120$	$3340 \pm180$	$0.15 \pm0.03$	$0.026 \pm0.02$	168	-616
Fe+Ni+Co		[Ni II]	$0.077 \pm0.009$	$-470 \pm940$	$3000 \pm1200$			(250)
Fe+Ni+Co		[Co III]	$0.97 \pm0.33$	$2300 \pm1500$	$2350 \pm660$
Fe+Ni+Co+Ca	[6800, 7600]	[Fe II]	$0.59 \pm0.01$	$680 \pm110$	$3360 \pm110$	$0.117 \pm0.024$	$0.035 \pm0.026$	171	-592
Fe+Ni+Co+Ca		[Ni II]	$0.069 \pm0.019$	$600 \pm280$	$4460 \pm830$			(252)
Fe+Ni+Co+Ca		[Co III]	$1.21 \pm0.35$	$980 \pm530$	$2820 \pm500$
Fe+Ni+Co+Ca		[Ca II]	$0.030 \pm0.019$	$-150 \pm590$	$2800 \pm1100$

We estimate the Ni-to-Fe abundance ratio following the approach of Maguire et al. [77]. Because Ni II and Fe II have similar ionization potentials, their number-density ratio is expected to approximately trace the elemental Ni/Fe abundance ratio [38]. Assuming optically thin emission, the ratio of the [Ni II] λ7378 and [Fe II] λ7155 line fluxes follows the equation:

where $L_{7378}$ and $L_{7155}$ are the measured flux ratio of [Ni II] 7378Å and [Fe II] 7155Å, respectively. $d_{\mathrm{cNiII}}/d_{\mathrm{cFeII}}$ is the departure coefficients from LTE for the appropriate states of Ni II and Fe II. Following Maguire et al. [77], we adopt a mean value $d_{\mathrm{cNiII}}/d_{\mathrm{cFeII}} = 1.8 \pm0.6$ estimated from the W7 model at +330d and allow the temperature to vary over 3000 − 8000K in nebular phase [38]. The range of the inferred $n_{\mathrm{NiII}}/n_{\mathrm{FeII}}$ for different line-identification models are listed in Table 3.

NIR transitions of iron-group ions can provide additional constraints on the intensity, velocity, and width of optical features arising from the same ions. Using a sample of nebular-phase optical+NIR X-Shooter spectra of SNe Ia, Flörs et al. [34] found that the flux ratio of the [Fe II] λ7155 and λ12570 lines follows an approximately log-linear temporal evolution with time since explosion:

with an intrinsic scatter of 0.06 dex. This trend may be explained by the natural cooling and dilution of the ejecta. Because of the relatively large uncertainties associated with the optical–NIR flux calibration and with the empirical relation itself, we do not impose this relation as a formal prior. Instead, we use it only as an independent consistency check on the inferred intensity and width of [Fe II] λ7155 in the 7200 Å line complex. From the flux calibrated FIRE+LRIS spectrum, we measure a line intensity ratio of $F_{12570}/F_{7155} = 0.76 \pm0.07$, in good agreement with the expected ratio between +350-380 days after explosion. In principle, the strength of Ni II λ7378 can similarly be constrained using Ni II λ19390 [34, 27]. In our case, however, the relatively low S/N and the strong telluric absorption blueward of Ni II λ19390 preclude a robust measurement of this feature.

Overall, the constant-continuum models are preferred, and the best-fit Fe+Ni+Co configuration yields a small $n_{\mathrm{NiII}}/n_{\mathrm{FeII}} = 0.026 \pm0.020$. Following the treatment in Flörs et al. [34], we account for the continued conversion of radioactive Co into Fe and scale the inferred abundance ratio at the observed epoch to its asymptotic value at $t \to\infty$. We thereby obtain an upper limit on the asymptotic Ni-to-Fe mass ratio of $M_{\mathrm{Ni}}/M_{\mathrm{Fe}} \lesssim0.05$. Adding Ca II to the line model is disfavored by the BIC and does not significantly improve the quality of the fit. Even when Ca II is included, the inferred abundance remains consistent with $M_{\mathrm{Ni}}/M_{\mathrm{Fe}} \lesssim0.06$ at $t \to\infty$. However, this result should be interpreted with caution. As shown in the comparison between different fitting settings, the best-fit result is highly dependent on the estimate of the continuum and the S/N because this line complex is highly blended. Nonetheless, the fittings we have largely exhausted the potential configurations and range of ion ratio, and we can conservatively conclude that the upper limit of $M_{\mathrm{Ni}}/M_{\mathrm{Fe}}$ is $\sim0.06$ at $t_{\infty}$.

both cases, the small Ni-to-Fe mass ratio is consistent with a sub-MCh Double Detonation progenitor model ([111, 34, 64]). The implications for the progenitor systems will be further discussed in § 5.

Discussion

Early-time observations of SN 2023bee revealed a pronounced flux excess within the first few days after explosion, with proposed origins including ejecta–companion interaction, circumstellar interaction, shallow or surface ${}^{56}\mathrm{Ni}$, and a sub-MCh double-detonation explosion ([51, 124]). Both Hosseinzadeh et al. [51] and Wang et al. [124] found that ejecta–companion interaction can provide a reasonable description of the early optical emission, with inferred binary separations of order $a \sim(1\text{--}3) \times10^{11}\,\mathrm{cm}$ under the assumption of a Roche-lobe-filling non-degenerate companion, corresponding approximately to a main-sequence companion rather than a red giant. However, Wang et al. [124] showed that the companion-interaction model does not simultaneously reproduce the full UV–optical evolution. They further explored alternative scenarios, finding that shallow ${}^{56}\mathrm{Ni}$ extending into the outer $\sim0.1\text{--}0.2\,M_{\odot}$ of the ejecta can reproduce the early optical rise, while a sub-MCh double-detonation model with a $\sim1.0\,M_{\odot}$ CO white dwarf and a thin $\sim0.01\text{--}0.02\,M_{\odot}$ He shell provides a comparably good description of the early light-curve morphology. Thus, although companion interaction remains a viable interpretation of the early excess, the early-time observations alone do not uniquely determine its physical origin and motivate independent constraints from the nebular phase.

In our late-time spectra, we find no significant narrow H or He emission associated with material stripped from a non-degenerate companion. Using the Botyanszki et al. [17] prescription, we derive nominal upper limits of $M_{\mathrm{stripped}}^{\mathrm{H}} \lesssim5.0 \times10^{-4}\,M_{\odot}$ and $M_{\mathrm{stripped}}^{\mathrm{He}} \lesssim1.2 \times10^{-4}\,M_{\odot}$. Accounting for the suppression of optical H emission by ejecta optical-depth effects following Dessart et al. [26] weakens the more conservative hydrogen constraint to $M_{\mathrm{stripped}}^{\mathrm{H}} \lesssim5.0 \times10^{-3}\,M_{\odot}$. These limits remain well below the stripped masses predicted by ejecta–companion interaction calculations. For example, Boehner et al. [14] found $\sim0.24\text{--}0.37\,M_{\odot}$ of stripped material for main-sequence companions with masses of $1.13\text{--}1.53\,M_{\odot}$ and binary separations of $2.13\text{--}6.58\,R_{\odot}$, partially overlapping the parameter space relevant to the companion-interaction interpretation of SN 2023bee. The observed limits therefore strongly disfavor such a classical single-degenerate configuration. Similar non-detections have been reported for several other SNe Ia with early flux excesses, including SNe 2017cbv, 2018oh, 2019yvq, and 2021aefx ([104, 29, 118, 112, 120, 50]). In light of these repeated non-detections in SNe Ia that otherwise show early-time flux excesses, companion interaction with a hydrogen- or helium-rich non-degenerate star is unlikely to be the dominant origin of early excesses in the population. More generally, systematic searches of nebular spectra have shown that narrow H/He emission is rarely detected in normal SNe Ia ([119]), disfavoring the SD progenitor scenario as the common origin of normal SNe Ia.

The nebular iron-group emission provides an independent constraint on the explosion conditions. At late phases, [Fe II] emission is dominated by the decay products of radioactive ${}^{56}\mathrm{Ni}$, whereas [Ni II] primarily traces stable Ni synthesized during the explosion. The relative strength of these species therefore probes the degree of neutronization and, consequently, the central density of the pre-explosion progenitor. Modeling of the blended $\sim7300$Å complex in SN 2023bee yields a weak [Ni II] contribution, with $L_{7378}/L_{7155} \sim0.1\text{--}0.2$ and a conservative asymptotic Ni-to-Fe mass-ratio limit of $M_{\mathrm{Ni}}/M_{\mathrm{Fe}} \lesssim0.06$. This places SN 2023bee toward the low-Ni end of the distribution measured in nebular samples ([77, 43]) and is more naturally reproduced by relatively low-density, sub-MCh burning than by standard high-density near-MCh delayed-detonation models ([34]). However, the inferred abundance ratio and physical interpretation remain subject to substantial systematic uncertainty. In particular, the method of Flörs et al. [34] does not account for variations in the chemical composition or density structure of the inner ejecta and assume that each ion has a free density. As Blondin et al. [10] pointed out, the intensity of [Ni II] is also dependent

on the ionization status and mixing of isotopes. For example, the lower density of [Ni II] in sub-$M_{\mathrm{Ch}}$ model suppresses the recombination process of Ni ions, resulting less Ni$^+$ ions and thus lower [Ni II] intensity compared to $M_{\mathrm{Ch}}$ explosion. In addition, the strong blending in the 7300Å region makes the result sensitive to the continuum placement, the adopted [Co III] and [Ca II] contributions, and the relatively modest S/N of the spectrum. We therefore regard the low Ni/Fe ratio primarily as evidence favoring lower-density burning rather than as a unique identification of a particular explosion mechanism.

Conclusion

We present late-time optical and NIR spectroscopy of SN 2023bee to investigate the progenitor and explosion scenario of this SN Ia with a prominent early flux excess. We find no significant narrow H or He emission associated with material stripped from a non-degenerate companion, placing upper limits of $M_{\mathrm{H}}^{\mathrm{stripped}} \lesssim5.0 \times10^{-3}\ M_\odot$ and $M_{\mathrm{He}}^{\mathrm{stripped}} \lesssim1.2 \times10^{-4}\ M_\odot$. These limits are difficult to reconcile with the non-degenerate companion configuration required by a simple ejecta–companion interpretation of the early excess. Modeling of the nebular $\sim7300$Å iron-group complex further yields a weak [Ni II] contribution and a conservative Ni-to-Fe mass-ratio limit of $M_{\mathrm{Ni}}/M_{\mathrm{Fe}} \lesssim0.06$, favoring a lower-density, sub-$M_{\mathrm{Ch}}$-like explosion, although the inference remains subject to systematic uncertainties in the modeling. Taken together, these late-time constraints provide a complementary view of the progenitor and explosion physics of SN 2023bee beyond that available from the early light curve alone. While no single diagnostic uniquely determines the progenitor scenario, the combined absence of stripped companion material and weak stable-Ni signature is consistent with a sub-$M_{\mathrm{Ch}}$ explosion and motivates consideration of a double-detonation origin for the early excess.

More broadly, nebular observations increasingly point toward a heterogeneous SN Ia population: stringent limits on stripped material place strong constraints on classical single-degenerate systems, while the diversity of nebular iron-group abundances indicates a range of central densities and explosion conditions. Recent nebular spectroscopy with JWST has revealed strong stable-Ni emission and signatures of high-density burning in several normal SNe Ia, including SNe 2021aefx, SN 2022aiq, and SN 2024gy, broadly consistent with near-$M_{\mathrm{Ch}}$ explosions [66, 3, 65] (DerKacy et al. 2023). In contrast, different analyses of SN 2022xkq have favored a near-$M_{\mathrm{Ch}}$ delayed-detonation explosion [24] and a sub-$M_{\mathrm{Ch}}$ explosion [65], respectively, which highlights the sensitivity of progenitor inferences from nebular spectra to the reduction method and spectroscopic modeling treatments. Such diversity suggests that multiple progenitor and explosion pathways likely contribute to the observed SN Ia population.

Establishing how these late-time properties relate to early light-curve behavior will require larger samples combining high-cadence observations during the first days after explosion with high-S/N nebular optical, NIR, and MIR spectroscopy. In particular, infrared observations provide access to less blended iron-group transitions and therefore cleaner constraints on stable isotopes and the physical conditions of the inner ejecta. Together with improved multidimensional non-LTE radiative-transfer calculations and atomic data, these observations will enable early-time signatures to be connected more directly with the composition, density, and geometry of the ejecta on an event-by-event and ultimately population-wide basis.

ACKNOWLEDGMENTS

This work was partially supported by NASA ADAP grant 80NSSC22K0494.

Q.W. is supported by the Sagol Weizmann-MIT Bridge Program. M.R.S. is supported by the STScI Postdoctoral Fellowship. P.C. acknowledges support from the Zhejiang Provincial Top-Level Research Support Program. Y.Z. acknowledges support from the MAOF grant (no. 12641898) and visitor support from the Observatories of the Carnegie Institution for Science, Pasadena, CA, where part of this work was completed. K.W.D. is supported in part by a University of California President’s Pre-Professoriate Fellowship (PPPF). The UCSC team is supported in part by NASA grants 80NSSC23K0301 and 80NSSC24K1411; and a fellowship from the David and Lucile Packard Foundation to R.J.F. A.G-Y.’s research is supported by ISF, IMOS and BSF grants, as well as the André Deloro Institute for Space and Optics Research, the Center for Experimental Physics, a WIS-MIT Sagol grant, the Norman E. Alexander Family M Foundation ULTRASAT Data Center Fund, and Magnus Konow in honor of his mother Olga Konow Rappaport; A.G-Y. is the incumbent of the Arlyn Imberman Professorial Chair. W.B.H. acknowledges support from the National Science Foundation Graduate Research Fellowship Program under Grant No. 2236415. C.D.K. gratefully acknowledges support from the NSF through AST-2432037, the HST Guest Observer Program through HST-SNAP-17070 and HST-GO-17706, and from JWST Archival Research through JWST-AR-6241 and JWST-AR-5441. C.L. is supported by DoE award # DE-SC0025599. G.N. is funded by NSF CAREER grant AST-2239364, supported in-part by funding from Charles Simonyi. G.N. also gratefully acknowledges NSF support from OAC-2311355, AST-2432428, as well as AST-2421845 and funding from the Simons Foundation for the NSF-Simons SkAI Institute. GN is also supported by the DOE through the Department of Physics at the University of Illinois, Urbana-Champaign. A.S. is is supported by a VILLUM FONDEN Villum Experiment grant (VIL69896) and a research grant (VIL54489) from VILLUM FONDEN.

The Young Supernova Experiment (YSE) and its research infrastructure is supported by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (ERC Grant Agreement 101002652, PI K. Mandel), the Heising-Simons Foundation (2018-0913, PI R. Foley; 2018-0911, PI R. Margutti), NASA (NNG17PX03C, PI R. Foley), NSF (AST-1720756, AST-1815935, PI R. Foley; AST-1909796, AST-1944985, PI R. Margutti), the David & Lucille Packard Foundation (PI R. Foley), VILLUM FONDEN (project 16599, PI J. Hjorth), and the Center for AstroPhysical Surveys (CAPS) at the National Center for Supercomputing Applications (NCSA) and the University of Illinois Urbana-Champaign.

Pan-STARRS is a project of the Institute for Astronomy of the University of Hawaii, and is supported by the NASA SSO Near Earth Observation Program under grants 80NSSC18K0971, NNX14AM74G, NNX12AR65G, NNX13AQ47G, NNX08AR22G, 80NSSC21K1572, and by the State of Hawaii. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, STScI, NASA under grant NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, NSF grant AST-1238877, the University of Maryland, Eötvös Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Some of the data presented herein were obtained at Keck Observatory, which is a private 501(c)3 non-profit organization operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

A major upgrade of the Kast spectrograph on the Shane 3 m telescope at Lick Observatory was made possible through generous gifts from the Heising-Simons Foundation as well as William and Marina Kast. Research at Lick Observatory is partially supported by a generous gift from Google.

This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.

This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NNA12R55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grant ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile.

Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grant Nos. AST-1440341, AST-2034437, and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, and IN2P3, France. Operations are conducted by COO, IPAC, and UW.

YSE-PZ was developed by the UC Santa Cruz Transients Team, supported in part by NASA grants NNG17PX03C, 80NSSC19K1386, and 80NSSC20K0953; NSF grants AST–1518052, AST–1815935, and AST–1911206; the Gordon & Betty Moore Foundation; the Heising-Simons Foundation; a fellowship from the David and Lucile Packard Foundation to R. J. Foley; Gordon and Betty Moore Foundation postdoctoral fellowships and a NASA Einstein fellowship, as administered through the NASA Hubble Fellowship program and grant HST-HF2-51462.001, to D. O. Jones; and a National Science Foundation Graduate Research Fellowship, administered through grant No. DGE-1339067, to D. A. Coulter.

Facilities: ATLAS, Keck:I, Magellan:Baade, Nickel, PO:1.2m, PS1, Shane, Swope

Software: astropy ([5, 6]), YSE-PZ ([57, 21] (Coulter et al. 2023)), Matplotlib ([53]), SciPy ([122]), NumPy ([44]), pyphot ([37]), emcee ([36]).