© The Authors 2024

1. Introduction

Supernovae (SNe), the luminous explosive deaths of stars, are an important ingredient in the chemical evolution of galaxies and an indicator of the final stages of the evolution of the exploding star. Core-collapse SNe (CCSNe), in particular, are the deaths of massive stars with zero-age main-sequence (ZAMS) masses of ≳8 M_⊙. CCSNe are divided into types Ib (H-poor), Ic (H- and He-poor), IIb (a small amount of H), and II (H-rich) based on the pre-SN structures and compositions of their progenitor stars (Filippenko 1997), which in turn is based on different kinds of mass loss during the lives of these stars. Mass loss is an especially murky aspect of stellar evolution (e.g., Smith 2014) that has a profound effect on the properties of the SN. SNe can thus be used to study the poorly understood late phases of the evolution of the massive progenitor stars.

Interaction with a binary companion can strip the hydrogen and/or helium envelope of a star, resulting in a stripped-envelope SN (SE-SN) of type IIb, Ib, or Ic depending on how much of the envelope is removed (Podsiadlowski et al. 1992; Yoon 2015; Drout et al. 2023). In addition to binary interaction, winds (at least those of very massive single stars; e.g., Yoon 2015) and pre-SN eruptions (e.g., Pastorello et al. 2007) can also strip the envelope and create circumstellar matter (CSM) which the ejecta of the SN can then interact with, resulting in extra luminosity and various signatures of interaction in light curves and spectra, such as narrow emission lines (e.g., in SNe Ibn and Icn, where “n” stands for narrow; Pastorello et al. 2008; Gal-Yam et al. 2022), but which mass-loss processes affect which progenitor stars is still somewhat unclear (Smith 2014).

Das et al. (2024) recently examined a sample of double-peaked SNe Ibc detected by the Zwicky Transient Facility (ZTF; Bellm et al. 2019; Graham et al. 2019). The first peak in these SNe is interpreted as the shock-cooling phase after a shock breakout from an extended envelope of material previously ejected in an eruption, with the second peak corresponding to the typical ⁵⁶Ni-powered peak of a SN Ibc. In superluminous SNe (SLSNe), early bumps in the light curve (e.g., Chen et al. 2023) have been suggested to result from such shock cooling (Piro 2015), the shock possibly driven by a magnetar (Kasen et al. 2016). Another double-peaked SN Ic was studied by Kuncarayakti et al. (2023), who attribute the first peak to ⁵⁶Ni decay and the second to circumstellar interaction (CSI) instead. Early peaks and bumps in SN light curves can also result from pre-SN eruptions, at least in H-rich interacting SNe IIn (e.g., in SN 2009ip and similar events; Mauerhan et al. 2013; Pastorello et al. 2018). In either case, double-peaked and other SNe with CSM-related features in their light curves can be used to probe the recent mass-loss history of the SN progenitor.

It is also possible that the later peaks in double-peaked events are not powered by CSI, but by a central engine. SN 2005bf (Anupama et al. 2005; Tominaga et al. 2005; Folatelli et al. 2006) is a double-peaked SN Ib that was suggested, among other things, to be powered by a double-peaked ⁵⁶Ni distribution (Tominaga et al. 2005) or a magnetar central engine (Maeda et al. 2007). Taddia et al. (2018) also suggest these scenarios, the latter possibly created by jets, to be responsible for the second peak in the SN Ic PTF11mnb (but the ⁵⁶Ni required an initial progenitor mass of 85 M_⊙, while the magnetar explanation could involve a less massive progenitor). SN 2019cad (Gutiérrez et al. 2021) is another double-peaked SN Ic in which a magnetar is considered a possibility for powering the second peak. In SN 2022jli, also of Type Ic, the favored power source is super-Eddington accretion from a companion star onto the newborn compact object instead (Chen et al. 2024). These examples demonstrate the diversity of double-peaked SE-SNe and their proposed mechanisms.

Here we study SN 2023aew, a remarkable seemingly type-changing double-peaked SE-SN with several peculiar properties. The discovery of SN 2023aew was reported by ZTF on 2023 Jan 23 (MJD = 59968; Munoz-Arancibia et al. 2023), and it was assigned the internal name ZTF23aaawbsc. However, even before this, on 2023 Jan 21 (MJD = 59965), it was also detected by the Gaia mission (Gaia Collaboration 2016) and assigned the Gaia Alerts¹ (Hodgkin et al. 2021) name Gaia23ate. The early evolution of the SN, including its initial peak, was missed by ZTF because of solar obscuration and because of technical problems and bad weather at Palomar Observatory (Sollerman, priv. comm.); the last Gaia non-detection is from 2022 Dec. 18 (MJD = 59932, i.e. over a month prior to discovery). However, after the original submission and arXiv release of this paper, Sharma et al. (2024) released another study of SN 2023aew, with some additional data including a serendipitous Transiting Exoplanet Survey Satellite (TESS) light curve that covers the first peak until shortly before the discovery.

The SN was initially classified as a SN IIb (Wise et al. 2023). However, after ∼90 d of slow decline, SN 2023aew exhibits a considerable rebrightening to a second peak brighter than the discovery magnitude, during which it resembles a SN Ic (Frohmaier et al. 2023; Hoogendam et al. 2023). The second peak is in itself both spectroscopically and photometrically peculiar for a SE-SN: the light-curve peak is broader than those of normal SNe Ic, even of “luminous” SNe (LSNe) of similar peak luminosity (Gomez et al. 2022), and the spectral evolution quite slow. The delay between the two peaks is, furthermore, very long compared to other double-peaked SE-SNe (e.g., Das et al. 2024).

In this paper, we present our photometric and spectroscopic follow-up observations of this remarkable transient in the optical and infrared. We analyze the light curve and spectral evolution, comparing both to other SNe Ic (especially double-peaked SNe), in order to assess numerous different physical scenarios. We present the observations we use, including proprietary and public ones, in Sect. 2, our analysis in Sect. 3, models of the photometry in Sect. 4, and a discussion of our results in Sect. 5. We present our conclusions in Sect. 6. Throughout this paper, magnitudes are in the AB system (Oke & Gunn 1983) unless otherwise noted, and we assume a ΛCDM cosmology with parameters H₀ = 69.6 km s⁻¹ Mpc⁻¹, Ω_M = 0.286 and Ω_Λ = 0.714 (Bennett et al. 2014). Reported uncertainties correspond to 68% (1σ) and upper limits to 3σ.

2. Observations and data reduction

We made use of the public gr-band ZTF light curve of SN 2023aew – through the Automatic Learning for the Rapid Classification of Events (ALeRCE) broker (Förster et al. 2021)² – and the public Gaia Alerts G-band light curve, and also retrieved the public classification spectrum of SN 2023aew (Wise et al. 2023) from the Transient Name Server (TNS³). We also extracted the public TESS light curve from Sector 60 and converted it to AB magnitudes (roughly equivalent to the i band) using TESSreduce⁴ (Ridden-Harper et al. 2021), which performs difference imaging, aperture photometry, and flux calibration based on the PANSTARRS1⁵ (Flewelling et al. 2020) catalog; the 200-s individual exposures were combined into 12-h bins. Finally, we obtained public near-ultraviolet (NUV) images from the Ultra-violet/Optical Telescope (UVOT) on the Neil Gehrels Swift Observatory (Gehrels et al. 2004) between MJD = 60076.5 and MJD = 60085.4 and a late-time template image from the NASA High Energy Astrophysics Science Archive Research Center (HEASARC) Data Archive⁶. The TESS and Swift data were first presented by Sharma et al. (2024). In addition to public data, we have obtained an extensive set of imaging and spectroscopic observations in the optical and infrared.

2.1. Imaging observations and photometry

The full logs of our photometric measurements, including the UVOT and TESS magnitudes, are shown in Tables A.1–A.10. Optical and near-infrared (NIR) imaging of SN 2023aew was done using the following instruments:

• the Alhambra Faint Object Spectrograph and Camera (ALFOSC) on the 2.56-m Nordic Optical Telescope (NOT) located at the Roque de los Muchachos Observatory, La Palma, as part of the NOT Unbiased Transient Survey 2 (NUTS2) collaboration⁷ (proposal 66-506, PIs Kankare, Stritzinger, Lundqvist): ugriz;

• the Las Cumbres Observatory Global Telescope (LCOGT) network (proposal OPTICON 22A/012, PI Stritzinger): uBgri;

• the 0.7-m GROWTH-India Telescope (GIT, Kumar et al. 2022) situated at Indian Astronomical Observatory (IAO), Hanle, India: griz. The data were obtained under regular proposals covering three cycles (C): GIT-2023-C1/C2/C3-P03 (PI: RS Teja);

• the IO:O instrument on the 2.0-m robotic Liverpool Telescope (LT) on La Palma (proposal PL23A18, PI Frohmaier): gri;

• the Near Infrared Camera Spectrometer (NICS; Baffa et al. 2001) on the 3.6-m Telescopio Nazionale Galileo (TNG) on La Palma (proposal A47TAC_37, PI Valerin): one epoch of JHKs;

• the Nordic Optical Telescope near-infrared Camera and spectrograph (NOTCam) on the NOT as part of NUTS2: JHKs;

• the Nishiharima Infrared Camera (NIC) on the 2.0-m Nayuta telescope at the Nishi-Harima Astronomical Observatory: JHKs.

The ALFOSC images were reduced using the custom pipeline Foscgui⁸, which performs bias subtraction and flat-field correction. The preprocessing and photometry of the GIT images were performed using a Python-based custom pipeline described in Kumar et al. (2022). LCOGT data were reduced by the automatic pipeline BANZAI⁹ which performs similar steps, also including astrometry calibration and bad-pixel masking. The LT images were similarly automatically reduced by the IO:O pipeline¹⁰. The Image Reduction and Analysis Facility (IRAF¹¹) package notcam¹² was used to perform sky subtraction, flat-field correction, distortion correction, bad-pixel masking, and co-addition of the NOTCam NIR images. TNG/NICS NIR images were reduced in a similar manner using IRAF tasks and the photometry performed using the SNOoPY¹³ pipeline. Similar reduction and point-spread-function photometry was performed for Nayuta/NIC data as well.

Aperture photometry was performed on the LCOGT, NOT, and LT images using Starlink¹⁴Gaia¹⁵ (Currie et al. 2014). The zero points of each image were calibrated using field stars from the Sloan Digital Sky Survey (SDSS) Data Release 17 (Abdurro’uf et al. 2022) in the ugriz filters, from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) in the JHKs bands, and from the AAVSO Photometric All-Sky Survey¹⁶ (APASS) data release 9 (Henden et al. 2016) in the B band. Color terms were taken into account for the LCOGT¹⁷ and LT¹⁸ photometry; zero color terms were assumed for NOT and ZTF photometry. Host template subtraction was not performed in the optical, as the location of the SN has very little contribution from the host. No source is detected at the location of the SN in pre-explosion PANSTARRS1 r-band images, down to a limiting magnitude of 23.2 mag.

For NUV photometry, we used the tasks uvotimsum and uvotsource in HEASOFT¹⁹ v.6.31.1 to co-add individual exposures and measure the source brightness, respectively. A 5-arcsecond circular aperture was used to extract the source flux, with background flux subtracted using a 10-arcsecond circular region; the same aperture was used in the NUV templates from 2023 Dec. 14 in order to subtract any host contamination from the earlier measured fluxes, similarly to Sharma et al. (2024).

In Fig. 1, we show a three-color (gri) image of the field around SN 2023aew, created using our NOT/ALFOSC images on MJD = 60134.0, when the seeing full width half-maximum (FWHM) was ∼0$\stackrel{″}{.}$4. This figure shows the location of the SN, seemingly in the outskirts of a spiral galaxy (SDSS J174049.60+661229.1) and close to a smaller galaxy (SDSS J174050.55+661220.7). These host galaxy candidates are examined in Sect. 3.1. We refer to the apparent nucleus of the spiral galaxy as Galaxy A, the smaller galaxy as Galaxy B, and the spiral galaxy itself as Galaxy C.

Fig. 1. Three-color (blue: g, green: r, red: i) cutout of the field of SN 2023aew, observed using NOT/ALFOSC on MJD = 60134.0. Nuclei of host galaxy candidates, Galaxy A and Galaxy B, are labeled. These galaxies are located at z = 0.30 behind the true host galaxy, Galaxy C (see Sect. 3.1).

2.2. Spectroscopic observations

We obtained a total of 26 optical spectra and one NIR spectrum of SN 2023aew during and after its second peak. The full log of spectroscopic observations is shown in Table A.11. These include

• eight spectra obtained using NOT/ALFOSC between MJD = 60061 and MJD = 60169 via NUTS2;

• six spectra between MJD = 60061 and MJD = 60113 using the fiber-fed KOOLS-IFU (Matsubayashi et al. 2019) integral field unit on the 3.8-m Seimei telescope (Kurita et al. 2020) at the Okayama Astrophysical Obsevatory in Japan (proposals 23A-K-0001, 23A-N-CT10, PI Maeda, and 23A-K-0017, PI Kawabata);

• four spectra between MJD = 60055 and MJD = 60075 with the SPectrograph for the Rapid Acquisition of Transients (SPRAT; Piascik et al. 2014) on the LT (proposal PL23A18, PI Frohmaier);

• four spectra between MJD = 60095 and MJD = 60126 with the Intermediate Dispersion Spectrograph (IDS) on the 2.5-m Isaac Newton Telescope (INT) on La Palma (proposal C54, PI Müller-Bravo);

• three spectra with the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) spectrograph on the 10.4-m Gran Telescopio Canarias (GTC), also on La Palma, between MJD = 60184 and MJD = 60227 (program GTCMULTIPLE2A-23A, PI Elias-Rosa);

• one spectrum on MJD = 60147 with the Gemini Multi-Object Spectrograph (GMOS) on the 8.1-m Gemini-North telescope at the Mauna Kea Observatories, Hawaii (program GN-2023A-Q-113, PI Ferrari);

• one NIR spectrum (11 700–24 700 Å) using TNG/NICS on MJD = 60064.

NOT/ALFOSC and GTC/OSIRIS spectra were reduced – bias-subtracted, flat-fielded, wavelength- and flux-calibrated, and corrected for telluric absorption – using Foscgui. LT/SPRAT spectra were automatically reduced by the LT pipeline (Barnsley et al. 2012). Performing the same steps, the TNG/NICS NIR spectrum was reduced through standard IRAF tasks. KOOLS-IFU spectra were bias-subtracted using IRAF, then reduced, sky-subtracted, and stacked using the Hydra package (Barden et al. 1994; Barden & Armandroff 1995). INT/IDS spectra were reduced using the package idsred²⁰. The GMOS spectrum was bias-subtracted, flat-fielded, and wavelength- and flux-calibrated using the gemini-gmos²¹ package in IRAF. All spectra of SN 2023aew we have obtained are available on the Weizmann Interactive Supernova Data Repository²² (WISeREP; Yaron & Gal-Yam 2012).

3. Analysis

3.1. Host galaxy and redshift

We used the SuperNova IDentification (SNID; Blondin & Tonry 2011) software to fit our spectra. The first spectrum was obtained with the low-resolution (R ∼ 100) Spectral Energy Distribution Machine (SEDM; Blagorodnova et al. 2018) on MJD = 59971.5, 18 rest-frame days after the first peak and 101 rest-frame days before the second (main) peak (see Sect. 3.2). It was reported by Wise et al. (2023) to resemble that of a SN IIb at z = 0.025. We find that this fit is only found by loosening the fit quality requirement parameter to rlapmin = 2.0 and forcing SNID to use type IIb, in which case the returned redshift is z = 0.026. However, the best matches from the GEneric cLAssification TOol (GELATO)²³ service using this redshift are also of type IIb, strengthening this association. By eye, the first spectrum of SN 2023aew 18 d after the first peak provides a reasonable match to SN 2003bg (Hamuy et al. 2009), a SN IIb, 20 d post-peak, as shown in Fig. 2. This SN is returned as the best fit by SNID. However, the SEDM spectrum has both a low resolution and a low signal-to-noise (S/N) ratio, making it not particularly constraining. The best SNID fits to our spectra taken during the second peak, on the other hand, are consistently those of post-peak SNe Ic²⁴ between z = 0.026 and z = 0.036. An exception is the first LT/SPRAT spectrum on MJD = 60055, which can fit either a post-peak SN Ib or a Ic in the same redshift range.

Fig. 2. Comparison of the SN 2023aew classification spectrum (red; Wise et al. 2023) and SN 2003bg (black; Hamuy et al. 2009), a SN IIb, at +20 d. The resemblance is fairly strong, but the SN 2023aew spectrum is quite noisy (uncertainties are shown as shaded regions). Telluric absorption lines (⊕) are marked with vertical gray lines.

The originally reported redshift of z = 0.025 indeed turns out to be correct. This conclusion is based on our GTC/OSIRIS spectrum from MJD = 60208.9, in which we find faint narrow (unresolved) emission lines that we identify as host contamination: Hα (z = 0.0252), Hβ (z = 0.0248), and [O III] λλ4959, 5007 (z = 0.0248 and 0.0247 respectively). Hα and [O III] λλ4959, 5007 are also seen in the MJD = 60226.9 spectrum, at z = 0.0254 and z = 0.0251 respectively. We also see lines at this redshift in our GTC/OSIRIS 2D spectra at the position of the large spiral galaxy close to the SN. We conclude that the spiral galaxy (Galaxy C in Fig. 1) is the host galaxy of SN 2023aew. In the following, we thus assume a redshift of z = 0.0250 ± 0.0003 for SN 2023aew, and all epochs, wavelengths, and absolute magnitudes are calculated with this assumption in mind. This redshift results in a luminosity distance of 109.8 ± 1.4 Mpc and a distance modulus of μ = 35.20 ± 0.03 mag.

Because of a lack of narrow emission lines from the host galaxy in our SN spectra before MJD = 60200 (see Sect. 3.3), we initially attempted to constrain the redshift of the SN by setting the position angle of the slit to intersect Galaxies A and B, which appeared to be the nuclei of the obvious host candidates. We positioned the slit on Galaxy A on MJD = 60068, 60077, and 60084 and on Galaxy B on MJD = 60117 and 60133. We then stacked the host galaxy spectra for total exposure times of 2700 s for Galaxy A and 4200 s for Galaxy B. The resulting combined spectra are shown in Fig. 3. It is evident that neither object has clear emission lines at the SNID-suggested redshift of SN 2023aew, z ≈ 0.03 (red dotted lines in Fig. 3). Instead, in both galaxies we do see clear emission lines of [O II] λ3727, [O III] λλ4959,5007, Hβ, and Hα – typical to star-forming SN host galaxies – but lines in both objects are at a redshift of z = 0.302 ± 0.001.

Fig. 3. NOT/ALFOSC spectra of Galaxies A (black) and B (blue) and the GTC/OSIRIS spectrum of SN 2023aew at +131 d showing narrow host lines (purple). Savitzky–Golay smoothing has been applied to the original ALFOSC spectra (gray). Identified lines of [OII], [OIII], and H at z = 0.30 are plotted as black dashed vertical lines. The same lines are also plotted at the redshift of SN 2023aew, z = 0.025, as red dotted vertical lines. No lines are detected at this redshift in Galaxies A and B.

The GTC/OSIRIS 2D spectrum from MJD = 60184.0 includes z = 0.30 emission lines from Galaxy A in addition to the lines from Galaxy C (while the spectrum at MJD = 60208.9 only includes z = 0.025 lines from Galaxy C, owing to a different slit position angle), and of the two sets of lines, the ones at z = 0.30 are considerably brighter. The fainter low-redshift lines from the true host galaxy are not seen in the earlier spectra (even stacked), most likely because of the lesser depth of the NOT/ALFOSC spectra. We thus infer that Galaxy A, which appears to be the nucleus of the host, Galaxy C, is, in fact, a separate galaxy shining through the disk of Galaxy C, most likely part of the same cluster as Galaxy B at z = 0.30.

3.2. Light curve

We performed numerical interpolation of the light curves using a Gaussian process (GP) regression algorithm (Rasmussen & Williams 2006) to obtain the peak magnitudes and epochs. We used the Python-based george package (Ambikasaran et al. 2015), which implements various different kernel functions. The widely used Matérn kernel with ν parameter of 3/2 was applied. The uncertainty in the peak epoch was estimated as the time range when the GP light curve is brighter than the 1σ lower bound on the peak brightness. The first peak is only covered by the TESS observations. We obtain a first-peak epoch of MJD${}_{\mathrm{TESS},\mathrm{peak}}=59953.4_{-2.4}^{+8.2}$. In the well-observed r band, we obtain a second-peak epoch of MJD_r, peak = 60074.9 ± 3.1. We use the latter as the reference date for all phases in this paper. We show the measured second-peak epochs and magnitudes in all bands in Table 1.

The ZTF discovery magnitude is r = 18.05 ± 0.10 mag, and the first g-band magnitude (∼2 d later) is g = 19.13 ± 0.12 mag, indicating a very red color soon after the first peak. The difference between the last TESS magnitudes, close to the i band in effective wavelength, and the first r points (r − TESS ∼ 0.5 mag, including a possible decline over 5 d) after the first peak, is also consistent with a very red color. We applied a Milky Way extinction correction of A_V, MW = 0.106 mag (Schlafly & Finkbeiner 2011), using the Cardelli et al. (1989) extinction law. The location of the SN at the outskirts of the host galaxy (see Fig. 1) and a lack of observed narrow Na I D absorption argue against a high extinction. Therefore we assumed a negligible host galaxy extinction, despite the red color after the first peak, which may be intrinsic to the event. This results in peak absolute magnitudes of M_TESS, peak = −17.88 ± 0.12 mag for the first peak and M_r, peak = −18.75 ± 0.04 mag for the more luminous second peak (the i band, more comparable to TESS data, peaks at a similar M_i, peak = −18.74 ± 0.04 mag). We note that Sharma et al. (2024) obtain different TESS magnitudes: their values are roughly 0.5 mag fainter than the ones returned by TESSreduce. Based on the g − r color of ∼1 mag, one might expect the TESS magnitudes to be brighter than in the r band, which is not the case in their Fig. 1, but such a difference can be seen in our values. The absolute-magnitude light curve, along with our GP fits, is shown in Fig. 4.

Fig. 4. NUV, optical, and NIR light curve of SN 2023aew. G-band errors are not reported at Gaia Alerts, and only the distance uncertainty is included. The shaded regions correspond to the 68% uncertainties of our GP fits. The decay rate of 56Co is indicated with a dashed black line. Our optical spectral coverage is indicated with gray vertical lines and our NIR spectrum with a black vertical line. The epoch of the classification spectrum (Wise et al. 2023) is marked with a blue vertical line. All magnitudes have been corrected for Galactic extinction. Phases refer to the second r-band peak (MJD = 60074.9).

The ∼90 d after the first peak show first a decline similar to the end of the ⁵⁶Ni bump and the beginning of the tail phase in a SE-SN, then a flattening inconsistent with ⁵⁶Co decay. The apparent undulations in the GP fit between −75 and −25 d are only a result of the slight over-fitting of the g-band light curve in this range (see Fig. 4): this is expected behavior, as there are few data points to anchor the fit and the evolution of the light curve in this range is slower than around and after peak (Stevance & Lee 2023). The rebrightening begins between −28 and −23 d, the last epoch of the plateau²⁵ and the first detection of the rebrightening, respectively. This is thus the constraint we obtain for the rise time to the second peak. As can be seen in Fig. 4 and in Table 1, in bands where we have enough data to reliably fit for the peak epoch – in the NIR and the NUV, there are not enough data points for this – the peak epochs are all consistent with the r band within 1σ. This is somewhat unusual if the second peak is a SN: in most SNe the bluer wavelengths tend to reach their peak first because of rapidly declining temperatures, but in SN 2023aew it is reached roughly simultaneously in all optical bands indicating an unusual temperature evolution.

After the second peak, the u and B bands decline linearly at 0.058 ± 0.002 and 0.043 ± 0.002 mag d⁻¹ (rest-frame) respectively. The NUV bands similarly decline quickly, and are faint and red (W1 − u ≈ 2 mag and W2 − W1 ≈ 1.5 mag) even near the peak. The g and especially the r and i bands initially decline more slowly, then show a steepening in the light curve around +40 d and a slight bump around +80 d. The bump is replicated in the J band as well. Such features could indicate the presence of CSM interaction or, since the gr-band bump corresponds roughly to the same brightness as the flattening after the first peak, possibly the continuing contribution of an earlier emission process. The post-peak decline is too swift for ⁵⁶Co decay, and a clearly ⁵⁶Co-powered tail phase is not reached during our observations. Our last photometry point at +257 d, consistent with the PANSTARRS1 detection limit, may be severely contaminated by host-galaxy emission (lines from the galaxy are detectable in our spectra ∼100 d before this; Sect. 3.1), and the decline rate at this point is likely underestimated. A faster decline, signifying dust formation or non-negligible gamma-ray leakage, is seen in, for example, SN 2020wnt (Gutiérrez et al. 2022; Tinyanont et al. 2023), but in SN 2020wnt the fully trapped decay is observed before this. Sometime between +110 and +135 d, the JHKs light curve starts rising again. Roughly simultaneously (between 130 and 140 d), the optical light curves briefly flatten. Even the NIR rebrightening has subsided by ∼190 d, at which point the decline has become faster in all bands. The gap in the NIR coverage prevents determining the duration of the rebrightening.

The color evolution of SN 2023aew, based on measured colors – when roughly simultaneous multi-band observations are available – and on our GP fits, is shown in Fig. 5. The g − r color, the only color available before the second peak, stays fairly constant at ∼1 mag from discovery until the sudden rise begins, at which point the SN becomes considerably bluer, reaching g − r ≈ 0.2 mag at peak. This rapid evolution to bluer colors and higher temperatures is reminiscent of the second peak of the double-peaked SN 2022xxf (Kuncarayakti et al. 2023). During the decline after the peak, the g − r color slowly reddens, reaching g − r ∼ 0.6 mag at +150 d. Meanwhile, the r − i color evolves even slower and becomes slightly bluer: at maximum light, r − i ≈ 0, but at +150 d, r − i ∼ −0.1 mag. This evolution is not typical for SE-SNe, which tend to quickly redden after the peak in all the optical colors (Stritzinger et al. 2018). The u − g color is more typical: it becomes bluer during the rise to the second peak, u − g ∼ 0.65 mag at its bluest, and subsequently reddens much faster than g − r, reaching u − g ∼ 1.6 mag around +30 d. The i − J and J − K colors are similar around the peak (∼ − 0.7 mag), but diverge around +40 d, with i − J reddening to ∼0.4 mag around +130 d and J − K only reaching ∼ − 0.3 mag.

Fig. 5. Measured colors of SN 2023aew in various bands when roughly simultaneous (Δt < 0.1 d) epochs are available (points) and subtractions between our GP fits of individual light curves when simultaneous photometry points at the different bands are scarce (lines and shaded regions).

We compare the light-curve behavior of SN 2023aew to other SE-SNe in Fig. 6. The first peak is similar in luminosity to SN 2003bg (Hamuy et al. 2009), but declines considerably slower from the beginning, both in and before the tail phase. The phase compared to SN2003bg also roughly matches that obtained from the SNID fit to the classification spectrum (see Fig. 2). The second peak, on the other hand, is compared to several SNe Ic. While the maximum luminosity of the second peak is close to that of SN 1998bw (Galama et al. 1998; McKenzie & Schaefer 1999), a broad-lined SN Ic (SN Ic-BL), and the LSNe²⁶ studied by Gomez et al. (2022), clear differences are apparent as well. During our observations, SN 2023aew does not reach the ⁵⁶Co-powered tail phase of SN 1998bw. The second peak of SN 2023aew evolves much slower than SNe Ic that attain a similar peak luminosity.

Fig. 6. Comparison of the light curve of SN 2023aew (points) to other SE-SNe (lines). The first peak at first resembles SN 2003bg (dashed purple line; Hamuy et al. 2009), until the flattening around −100 d. The second peak has a luminosity similar to the Type Ic-BL SN 1998bw (Galama et al. 1998; McKenzie & Schaefer 1999) and LSNe (Gomez et al. 2022) (solid lines), but a much slower evolution. Of the double-peaked SE-SNe (DP; dotted lines) pictured here – SNe 2019cad (Gutiérrez et al. 2021), 2020bvc (Ho et al. 2020; Izzo et al. 2020), 2022xxf (Kuncarayakti et al. 2023), and 2022jli (Chen et al. 2024) – only SN 2020bvc resembles the late-time light curve of SN 2023aew, and the separation between peaks is much shorter in all of them than in SN 2023aew.

We also show a comparison to four other double-peaked SNe Ic: SNe 2019cad (Gutiérrez et al. 2021), 2020bvc (Ho et al. 2020; Izzo et al. 2020), 2022jli (Chen et al. 2024), and 2022xxf (Kuncarayakti et al. 2023). All show a much shorter delay and/or a smaller contrast between the peaks than SN 2023aew. SNe 2019cad and 2022xxf evolve much faster than the second peak of SN 2023aew, while SN 2022jli shows an even slower, linear decline and periodic undulations that are not seen in SN 2023aew (an ∼80-day period is technically possible, but with only two bumps we cannot say this for certain). All three of these double-peaked SNe Ic have a different light-curve evolution and different proposed physical scenarios, but none of them resemble SN 2023aew. SN 2020bvc does show a similar peak luminosity and bumpy light curve as SN 2023aew, but its light-curve peak is again narrower and its spectrum that of a SN Ic-BL, with hints of an accompanying long gamma-ray burst that remained undetected (Ho et al. 2020; Izzo et al. 2020).

The double-peaked SE-SNe studied by Das et al. (2024) mostly do not closely resemble the light-curve evolution of SN 2023aew either; a similar delay between peaks is only seen in one event in the ZTF sample, and in this case (SN 2021uvy) the rise to both peaks is much slower and the “dip” between the peaks much less pronounced. However, the objects in Das et al. (2024) do exhibit a wide range of decline rates.

3.3. Spectroscopy

We present our optical spectra of SN 2023aew and the Wise et al. (2023) SEDM spectrum of the first peak in Fig. 7. Here we describe the main features of the spectrum and their evolution.

Fig. 7. Spectral evolution of the first (blue; Wise et al. 2023) and second peak of SN 2023aew. Different instruments are color-coded, and epochs refer to the second peak. Telluric absorption lines (⊕) are marked with solid light gray lines. Dash-dotted lines are at zero velocity (of the stronger line in case of a close doublet); long dashes at −2500 km s−1; short dashes at −6500 km s−1; and dotted lines at −8500 km s−1. When Savitzky–Golay smoothing has been applied, the original spectrum is plotted in gray.

A clear change is apparent between spectra taken during the first and second peak. The classification spectrum at −101 d relative to the second peak is characterized by lines of He I, Ca II and Hα (possibly contaminated by Si IIλλ6347,6371). The He Iλ5876 line is seen in absorption around −9000 km s⁻¹, and the Ca II NIR triplet absorption exhibits a similar velocity, while the 6678 and 7065 Å lines are dominated by emission. These features are similar to the Type IIb SN 2003bg around 10–20 d past maximum light (Hamuy et al. 2009; Mazzali et al. 2009), which is also suggested by the SNID fit at this epoch. The absorption minimum of the Hα line is at −12 500 km s⁻¹, assuming no contribution from Si IIλλ6347,6371 – a line dominated by Si II is unlikely, as its velocity would be only ∼−3000 km s⁻¹. A comparison of the pseudo-equivalent width of the Hα+Si IIλλ6347,6371 feature (∼160 Å) to the ones presented in Holmbo et al. (2023) indicates a type IIb and an Hα-dominated feature as well, while the strength of the He Iλ5876+Na Iλλ5890,5986 feature is more typical of SNe Ic. We note, though, that the spectrum of SN 2003bg (see Fig. 2) differs from prototypical SNe IIb (e.g., Filippenko et al. 1994; Ergon et al. 2014) in that the Hα line is still stronger than He lines at this phase. Some C IIλ6580 contribution cannot be excluded, though.

In the −19 d spectrum, taken only a few days after the start of the rebrightening, and other spectra taken before and at the second peak, we see a different set of lines, as indicated by the reasonable SNID matches to SNe Ic (Sect. 3.1)²⁷. An absorption feature close to Hα and an emission feature close to He Iλ7065 are still present, but the He Iλ6678 line (unfortunately disturbed by a telluric line) seems weak or nonexistent and He Iλ4471 is not visible. Lines of O I and Mg II, and the Ca II H&K lines, which are not (clearly) present in the −101 d spectrum, also appear, as does [Ca II] λλ7291,7323 possibly blended with [O II] λλ7320,7330. The blue part of the spectrum is now dominated by a plethora of iron lines (Dessart et al. 2021, 2023). Some of the features visible at −14 d could be lost in the low S/N of the −101 d spectrum, but the disappearance of He Iλ6678 and the appearance of O I lines do not seem to be attributable to this.

The He Iλ5876 absorption is now blended with the Na Iλλ5890,5986 doublet, and the absorption trough likely has a contribution from both. The velocity of the absorption minimum is ∼ − 6800 km s⁻¹ if dominated by Na I, or ∼ − 6000 km s⁻¹ if dominated by He I. This velocity is comparable to the absorption minimum of the Ca II NIR triplet at −14 d, ∼ − 6000 km s⁻¹, and that of the strong O Iλ7774 line, ∼ − 5500 km s⁻¹. If one assumes the absorption line close to Hα is still Hα, its velocity is ∼ − 11 000 km s⁻¹ (a slightly higher velocity of ∼ − 11 500 km s⁻¹ is obtained assuming C IIλ6580). This absorption and the accompanying emission peak might be a combination of Hα and the Si IIλλ6347,6371 feature, but Hα emission would then be much more blueshifted than other emission lines, suggesting a different origin. Spectral models by Dessart (2019) and Dessart et al. (2023) show blends of Fe II emission lines at 6148, 6149, 6238, 6247, 6417, 6433, 6456, and 6516 Å. These could produce the observed features, which peak at 6190 and 6480 Å.

Spectroscopic evolution is, afterwards, quite slow, and the velocities of the Na/He and O I lines decline only by about 1500 km s⁻¹ from −19 d to +148 d. The slow evolution is emphasized by the fact that the spectra of the second peak consistently resemble SNe Ic at post-peak epochs, even before this peak. The main changes between these epochs include the appearance of low-velocity (−2500 km s⁻¹) Sc II absorption lines; the gradual appearance of strong emission features at Na Iλλ5890,5986, [O I] λλ6300,6364, and at ∼6550 Å; and the fading of Mg IIλ7887 while Mg IIλ8224 persists. The initially iron-dominated peak at ∼4600 Å gradually moves slightly bluewards as it starts being dominated by Mg I] λ4571. Mg I lines around 5170 Å also appear in a P Cygni feature around the same time, with an absorption minimum at ∼ − 3000 km s⁻¹.

The ∼6550 Å feature appears roughly simultaneously with [O I] λλ6300,6364 and Mg I] λ4571 while the Fe II features in these regions fade. An obvious possibility for this line would be Hα; by definition this line is not seen in a normal SN Ic, but could be a result of interaction. At such a late stage of evolution, C IIλ6580 is unlikely. Alternatively, if the first peak is a true SN IIb (see Sect. 5), the [O I], Mg I], and ∼6550 Å features appear ≲200 d after that explosion, when forbidden/nebular lines from the first SN could be seen. The timing also coincides with the bump phase where the SN reaches the same brightness as the flattening between the two peaks (Sect. 3.2) and where nebular features from the first SN might start contributing noticeably if the flat underlying light curve has continued.

Such a strong late-time Hα line is not typical for a SN IIb, and [N II] is in any case considered the more likely origin by Jerkstrand et al. (2015) and Fang & Maeda (2018). Hα from CSI was seen in late-time spectra of the prototypical SN 1993J, but it was also much weaker than the line seen in SN 2023aew, and the shape of its profile was boxy as expected from a CSM shell (e.g., Patat et al. 1995). [N II] is typically very weak in late-time SN Ic spectra compared to [O I], but can be comparable to [O I] in SNe IIb with low progenitor masses. A [N II]/[O I] ratio of ∼2/3 is seen in SN 2023aew at +148 d; a 12-M_⊙ progenitor star model can produce such a ratio (Jerkstrand et al. 2015), and Fang et al. (2019) report some SNe IIb with similar ratios, though such SNe also tend to exhibit [O I]/[Ca II] ratios below 1, whereas in SN 2023aew this ratio is ∼1.7 at +148 d. [N II] is produced in the helium layer of the ejecta, where CNO burning results in a high nitrogen abundance (Jerkstrand et al. 2015), which is compatible with the helium contribution in the relatively early spectra. Since the line is quite strong for [N II] given the high [O I]/[Ca II] ratio, a combination of [N II] λλ6548,6584 and Hα is the most likely option; similar cases of Hα “contamination” are noted by Fang et al. (2019).

Fig. 8 shows the evolution of the Na/He absorption, the putative [N II] and [O I] lines, and the [Ca II] λλ7291,7323 + [O II] λλ7320,7330 feature. At late times (>  + 130 d), all of these features evolve to show a slightly double-peaked profile (i.e., two overlapping doublets at different velocities). The widths of the three features are similar (∼5000 km s⁻¹ after +130 d), and they seem to exhibit additional peaks blueshifted by ∼1500 km s⁻¹ relative to each line. The locations of these peaks compared to the other lines are consistent with an [N II] origin for the ∼6550 Å feature. Furthermore, if Fe II lines were still responsible for the putative [O I] and [N II], this would result in a redshift of several thousand km s⁻¹ for the peaks of these blends, which is not seen in other spectral lines. We thus consider it more likely that these emission features are correctly identified as [O I] λλ6300,6364 and [N II] λλ6548,6584 (likely with Hα contamination). We do note, though, that the latter is difficult to verify as the [N II] λ5754 feature falls on the absorption trough of the Na/He absorption and is not clearly seen.

Fig. 8. Left: evolution of the Na Iλλ5890,5986 line (dashed gray lines), possibly blended with He Iλ5876. The approximate evolution of the absorption minimum from a Gaussian fit after −19 d is marked with a solid gray line. The velocity evolution is slow during the second peak and the subsequent decline. Middle: evolution of the putative [O I] λλ6300,6364 (dashed red lines) and [N II] λλ6548,6584 (dashed gray lines) doublets after the second peak. Both start to appear between +41 and +71 d and strengthen over time. Right: evolution of the [Ca II] λλ7291,7323 (dashed gray lines) + [O II] λλ7320,7330 (dashed red lines) feature. Zero velocity refers to the stronger line of each doublet. All forbidden features shown here are double-peaked at late times.

We present our NIR spectrum of SN 2023aew at −10 d in Fig. 9. This spectrum is rather noisy, but we see some features of the same elements as in the optical. If the absorption feature at ∼12 200 Å is identified with Paβ, its velocity is ∼ − 12 500 km s⁻¹, close to the putative Hα velocity and suggesting a hydrogen origin for that line, but other NIR lines of hydrogen are not clearly seen – Paα is located in a gap in the spectrum, but Brγ may be lost in the noise. Absorption consistent with He Iλ20 581 at ∼ − 8500 km s⁻¹ is seen in the spectrum as well as the emission of Mg IIλ21 369. We also see a possible C Iλ16 890 emission line, perhaps accompanied by weak P Cygni absorption. Both putative emission lines are blueshifted by ∼1500 km s⁻¹. These lines are not seen in the NIR spectra of typical SNe Ic shown by Shahbandeh et al. (2022), which instead show Mg Iλ14 878 and C Iλ14 543. SNe IIb do show Paβ and He Iλ20 581, but the latter mostly in absorption. As the spectrum is quite noisy, these identifications are uncertain.

Fig. 9. TNG/NICS NIR spectrum of SN 2023aew at −10 d. Possible features of H, He, Mg and C are seen, but the spectrum is noisy and atypical for SE-SNe. Line velocities are as in Fig. 7, with the addition of −12 500 km s−1 as a solid line. Savitzky–Golay smoothing has been applied.

For comparison between the second peak of SN 2023aew and normal SNe Ic, we obtained publicly available spectra of SNe 1994I (Filippenko et al. 1995; Modjaz et al. 2014) and 2007gr (Valenti et al. 2008; Modjaz et al. 2014; Shivvers et al. 2019) from WISeREP. We show this comparison at different epochs in Fig. 10. It is clear from this figure that SN 2023aew spectroscopically evolves slower than either comparison SN. Even the first spectrum of SN 2023aew, 19 d before the second peak, resembles that of SN 2007gr at much later times, with much stronger features. After the second peak, and especially after a few weeks post-peak, SN 2023aew begins to spectroscopically resemble SN 2007gr, with similar strengths and widths in the absorption lines; however, the broad emission peak at ∼6500 Å is not seen in late-time (> 60 d) spectra of either SN 1994I or SN 2007gr. Instead, only a narrow Hα line from the host galaxy is seen in the latter two. The resemblance to SN 2007gr is generally stronger than SN 1994I in terms of velocity and line strength. SN 2007gr still tends toward higher velocities than SN 2023aew, while SN 1994I exhibits even higher velocities with fewer and, until late times, generally weaker lines. The velocity comparison is illustrated in Fig. 11, where we show the evolution of the Na Iλλ5890,5896 absorption minimum.

Fig. 10. Spectra of SN 2023aew compared to well-observed ordinary SNe Ic. While neither SN 1994I (thick blue line) nor SN 2007gr (dashed red line) matches the early-time features of SN 2023aew, the latter does resemble SN 2007gr at later epochs. The feature at ∼6500 Å is not seen in late-time spectra of either SN 1994I or SN 2007gr. Overall the evolution of SN 2023aew is relatively slow.

Fig. 11. Evolution of the Na Iλλ5890,5896 absorption minimum velocity in SN 2023aew (black) after the rebrightening, compared to SNe 1994I (blue) and 2007gr (red). Uncertainties are estimated at 5%. For SN 2023aew, the phase refers to the second peak.

This is somewhat similar to what was seen in the unusual SLSN I, SN 2020wnt (Gutiérrez et al. 2022), which also showed a similarity to SN 2007gr, but with slower evolution and a ∼6500 Å feature. This object also shared some light-curve features in common with the second peak of SN 2023aew (Sect. 3.2). The ∼6500 Å feature was attributed to [Fe II] λ6456 and λ6518 in SN 2020wnt, but this is unlikely in SN 2023aew, as the line profile would have to be redshifted by ∼4000 km s⁻¹.

We have also compared SN 2023aew to other double-peaked SE-SNe: SN 2019cad (Gutiérrez et al. 2021), SN 2022jli (Chen et al. 2024; Moore et al. 2023), and SN 2022xxf (Kuncarayakti et al. 2023). We show this comparison in Fig. 12, where all epochs refer to the second peak, which is the brighter one in all objects. These comparison objects, as stated in Sect. 3.2, are fainter than SN 2023aew, except SN 2022xxf, and exhibit a shorter delay between the two peaks. In general, once again, we see a faster evolution in the spectra as well. However, some similarities are also seen: the last spectrum of SN 2019cad, especially, where the 6500 Å feature also appears (albeit weaker and unidentified), matches the spectrum of SN 2023aew well.

Fig. 12. Spectra of SN 2023aew (black) compared to double-peaked SE-SNe SN 2019cad (purple), SN 2022jli (thick green line) and SN 2022xxf (dashed brown line). All epochs refer to the second peak. SN 2023aew evolves slower than SNe 2019cad or 2022xxf and exhibits lower velocities, but an excellent match is seen with the latest spectrum of SN 2019cad. The 6500 Å feature, in particular, is replicated. The match with SN 2022xxf similarly improves at late times, although this feature does not appear. SN 2022jli provides a good match at ∼30–40 days and evolves slower than other comparison objects, but lacks a strong [O I] doublet at late times.

4. Modeling

4.1. Blackbody fits

We used the GP-interpolated light curves (see Sect. 3.2) to construct spectral energy distributions (SEDs) at roughly 10-day intervals between the first JHKs-band epoch (−11.4 d) and the last JHKs epoch at +135.6 d. We note, however, that we do not have NIR data at the optical peak, and close to that epoch, the NIR fluxes are most likely underestimated by the GP fit, leading to overestimated temperatures. Conversely, in the epochs after +70 d we lack the uB bands, and temperatures can be affected from the other direction. Additionally, we used the preceding epochs where only ZTF gr-band data are available and where the parameters are thus highly uncertain. We show the blackbody function at each epoch compared to the interpolated SED in Fig. 13. Line blanketing in the NUV results in clear deviation from a blackbody, and the NUV bands were ignored in the fits.

Fig. 13. Blackbody fits (lines) to the observed SEDs (points) at roughly 10-day intervals. NUV points (light gray) at +2 and +10 d clearly deviate from a blackbody because of line blanketing, and were ignored in the fits.

4.1.1. Evolution of the blackbody

The evolution of the blackbody parameters, the temperature and radius, is shown in Fig. 14. Initially, in the epochs where we can only use gr photometry, the color stays constant implying a temperature around 4000 K, while the radius decreases from 3 × 10¹⁵ cm to about 2 × 10¹⁵ cm close to the rebrightening. In the well-measured epochs, the radius stays fairly constant around 2 × 10¹⁵ cm, until a decline around +40 d, which coincides with the steepening of the ri-band light curve (Fig. 4). The radius then appears to slowly decline from 2 × 10¹⁵ cm at +40 d to ∼6 × 10¹⁴ cm at +120 d. Meanwhile, the temperature seems to reach a peak of ∼8100 K at the optical light curve peak (but note that this is where the NIR interpolation is likely increasing the temperature) before declining slowly to ∼5500 K at +40 d, after which it stays roughly constant. The last epoch at +136 d seems to exhibit a lower temperature, around 4700 K, and an expanding radius, but at late times the SED deviates somewhat from a blackbody, and the resulting parameters may be unreliable.

Fig. 14. Evolution of the blackbody parameters over time. The temperature from the fit rises during the rebrightening from a preceding plateau and is highest at the peak – albeit likely slightly overestimated – then declines to ∼5500 K and stays constant until ≳120 d. The radius slowly declines almost throughout the observed evolution.

The mostly continuous decrease of the radius seems to indicate that the photosphere moves inward through our observations. The rebrightening is mostly caused by an increase in temperature, with possibly a pause in the shrinking of the photosphere as well; hence, the blue bands did not peak before the red ones (Sect. 3.2). Even if one ignores the fits using only gr data, a radius of ∼2 × 10¹⁵ cm is already reached at −11 d. The rebrightening starts between −28 d and −23 d; if this corresponds to the explosion date of a SN, the ejecta velocity required to reach the blackbody radius ranges from ∼14 000 to ∼19 000 km s⁻¹. Such velocities are not supported by our spectra, where the absorption minimum of the Na Iλλ5890,5896 doublet is at 6800 km s⁻¹ at −19 d (or, if contaminated by He Iλ5876, even lower). If the feature at ∼6500 Å in the −101 d spectrum is Hα as the IIb classification suggests, the velocity of the absorption minimum is −12 500 km s⁻¹. As photospheric velocities tend to be smaller than the Hα velocity, reaching a radius of ∼3 × 10¹⁵ cm would take ≳28 d, which is compatible with the epoch of the spectrum ∼34 d after the explosion leading to the first peak.

4.1.2. Bolometric light curve

We used SuperBol (Nicholl 2018) to construct the bolometric light curve of the second peak of SN 2023aew, where multi-band photometry is available. The r band was used as the baseline, and the optical and NIR light curves in other filters were interpolated to the r-band epochs using a polynomial fit within SuperBol. Extrapolation of the very early and/or late light curve was performed by assuming a constant color, as polynomial fits often diverged wildly from the optical light curves. NUV data were not included due to the short time baseline that would necessitate extrapolation over most of the light curve, but SuperBol blackbody fits include simulated line blanketing below 3000 Å to account for the observed deviation of the NUV fluxes from a blackbody. We constructed both a pseudo-bolometric light curve using only the observed bands and a bolometric light curve corrected for missed emission using SuperBol’s blackbody fit. Additionally, we used our gr-band blackbody fits to estimate the bolometric luminosity before the rebrightening. These are displayed in Fig. 15, along with comparison events. Out of these, only SN 1998bw is more energetic than the second peak because of its slower decline.

Fig. 15. SuperBol (Nicholl 2018) pseudo-bolometric light curve of the second peak integrated over observed bands (black open circles), and a bolometric light curve based on blackbody fits (red filled circles), highly uncertain before the rebrightening. The difference is notable only around maximum light. The late-time decline is inconsistent with 56Co decay (solid black line). Pseudo-bolometric light curves of SNe 1998bw (Patat et al. 2001), 2019cad (Gutiérrez et al. 2021), 2022jli (Moore et al. 2023), and 2022xxf (Kuncarayakti et al. 2023) are included for comparison.

The blackbody correction is small except near maximum light, when the temperature also peaked (see Sect. 4.1.1). The evolution of the bolometric light curve is very similar to that of the r-band light curve, with bumps around 80 and 140 days – this is expected, as the bumps are replicated in the NIR as well. A tail phase with a decline consistent with the decay of ⁵⁶Co at full gamma-ray trapping is not reached during our observations. We also integrated over both light curves to determine the total radiated energy during the second peak. In the observed filters, we obtain a total of (3.67 ± 0.01)×10⁴⁹ erg, while the blackbody-corrected total radiated energy is (4.34 ± 0.04)×10⁴⁹ erg.

4.2. Light-curve fitting

We used the publicly available code Modular Open Source Fitter for Transients (MOSFiT²⁸; Guillochon et al. 2018) to fit the second peak of the light curve, assuming it was a H-poor SN. The code includes the ⁵⁶Ni decay model (Arnett 1982; Nadyozhin 1994), labeled default; a CSI model (Chatzopoulos et al. 2013; Villar et al. 2017; Jiang et al. 2020), labeled csm; a fallback accretion model (Moriya et al. 2018a), labeled fallback; and a magnetar central engine model with line blanketing correction (Nicholl et al. 2017), labeled slsn. These models are fitted using the dynamic nested sampling package dynesty.²⁹ We ran each fitting process until convergence, which typically took between 25 000 and 30 000 iterations.

We set simple uninformative uniform or log-uniform priors for each free parameter, summarized in Table 2 for each model. We set the lower limit on the characteristic ejecta velocity to be close to the Na I absorption minimum in our earliest spectrum, and the upper limit at twice the lower limit. An upper limit for the column density of neutral hydrogen in the host galaxy, n_H, host, was set at 2 × 10²¹ cm⁻² (corresponding to A_V ≲ 1 mag according to Güver & Özel 2009), as the host extinction is likely to be low. Other parameters common to all models include the explosion time before observations t_expl, the minimum temperature T_min, and the ejecta mass M_ej. The nickel model additionally includes the nickel fraction in the ejecta f_Ni, the opacity κ, and the opacity to γ rays κ_γ. The magnetar model, on the other hand, includes κ, the spin period P_spin, the magnetic field perpendicular to the spin axis B_⊥, the neutron star mass M_NS, and the angle between the magnetic field and spin axis θ_PB. In the CSI model we include the parameter s (where the CSM density as a function of distance behaves as ρ ∝ r^−s), which we fix to s = 0, indicating a CSM shell, and s = 2, indicating a wind. In addition, this model includes the minimum inner radius of the CSM, R₀, the CSM mass M_CSM, and the CSM density at R₀, ρ. We fix the density power-law parameters in the inner (ρ_ej ∝ r^−δ) and outer (ρ_ej ∝ r⁻ⁿ) ejecta, δ = 1 and n = 7, respectively (corresponding to normal values in H-poor ejecta: Chatzopoulos et al. 2013; Chevalier 1982). Finally, the fallback model includes κ, the linear-to-power-law accretion transition time t_tr, and the accretion luminosity at the transition time L₁. In total, the default, fallback, and csm models have 9 free parameters and the slsn model has 12. These numbers include a nuisance parameter σ, which describes the added variance required to match the model being fitted.

We also used MOSFiT to fit the light curve before rebrightening, ignoring the flattening after MJD = 60000, in order to obtain an estimate of the parameters of the first peak, again with the assumption that it was a genuine SN. In this case, we used only the default model, with the same priors as for the second peak, except the ejecta velocity prior was set as [8000:15 000] km s⁻¹, as informed by the −101 d spectrum. TESS data points were approximated as i-band photometry in this fit.

The results (median parameters and 1σ uncertainties) of our fits are listed in Table 3, and the model light curves are shown in Fig. 16. The corner plots of the fits are shown in Figs. B.1–B.6. All five models have trouble reproducing the second peak. Specifically, the main problem is that all optical bands show a simultaneous light-curve peak (Sect. 3.2), while in the MOSFiT models the blue bands, especially u, peak first. The u-band light curve shape is also difficult to reproduce. In the NUV, all our best fits include a much brighter peak prior to the observed epochs, which probably also influences the optical peaks. The ejecta nickel fraction f_Ni in the ⁵⁶Ni decay model is practically 1, which further disfavors this model for the second peak. In addition, the fallback model results in a t^−5/3 decline, which overpredicts late-time fluxes, disfavoring this model as well. The magnetar model provides the best fit by eye and has the highest likelihood score returned by MOSFiT, corresponding to the logarithm of the Bayesian evidence Z, favoring it over the other models, but we note that the κ and κ_γ parameters of this model cluster at the edges of their priors.

Fig. 16. MOSFiT light curve fits to the first and second peak of SN 2023aew using the different models described in the text. The flat part of the light curve after MJD = 60000 was ignored in the first-peak fit, as the 56Ni model does not include such features. All models have trouble reproducing the simultaneity of the second peak in all optical bands, specifically failing to reproduce the u-band peak, but by eye and by MOSFiT score, the magnetar model provides the best fit.

The ⁵⁶Ni model fit to the first peak is fairly good; unlike the second peak, the features within the fitted time range are reproduced by the model. The median ejecta and Ni masses are ∼3.3 and ∼0.2 M_⊙, respectively – fairly close to those found for SN 2003bg (∼4 and ∼0.2 M_⊙, respectively, according to Mazzali et al. 2009), which also exhibits a similar spectrum and peak luminosity. The caveat here is, of course, that the following flattening is not typical for a SN IIb.

5. Discussion

We now discuss the features of SN 2023aew in the context of other double-peaked SE-SNe, and several possible scenarios that might conceivably explain these features. These scenarios include two separate SNe – whether related to each other or not – that are each responsible for one of the peaks and happen on the same line of sight; one SN preceded by an eruption; one SN followed by delayed CSI; one SN followed by delayed activity of a magnetar central engine; and one SN followed by accretion onto a newborn compact object.

5.1. SN 2023aew vs. other double-peaked SE-SNe

Among the double-peaked SNe we have compared SN 2023aew against, it stands unique (though we note that at least one similar object, SN 2023plg, exists according to Sharma et al. 2024, but no spectra of its first peak exist). Other double-peaked SE-SNe, such as SNe 2019cad (Gutiérrez et al. 2021), 2020bvc (Ho et al. 2020; Izzo et al. 2020), 2022jli (Moore et al. 2023; Chen et al. 2024), and 2022xxf (Kuncarayakti et al. 2023), and the ZTF sample of double-peaked SNe Ibc from Das et al. (2024), all show a shorter separation between the light-curve peaks and/or a much less pronounced second peak. Apart from SN 2022jli, they also show a faster spectroscopic and photometric evolution than SN 2023aew, as do normal single-peaked SNe Ic. Out of these, the SN whose light curve most resembles SN 2023aew is SN 2020bvc – the peak is less broad, but the peak luminosity and late-time bumpy light curve are similar – but it shows much higher velocities and is classified as a SN Ic-BL.

Spectroscopically the best match to the late-time spectra of SN 2023aew is found in the Type Ic SN 2019cad: all features including the possible [N II] feature³⁰ at ∼6550 Å are replicated. However, this is only true in one epoch (46 d after the second peak or 88 d after explosion; Gutiérrez et al. 2021), before which SN 2019cad evolves faster than SN 2023aew. Later spectra of SN 2019cad are, unfortunately, not available. The ∼6550 Å feature is not seen in the normal SNe Ic we have examined (Filippenko et al. 1995; Hunter et al. 2009; Modjaz et al. 2014), nor in SN 2003bg at late times (Hamuy et al. 2009), nor in the double-peaked SN 2022xxf. SN 2022jli does show a feature at ∼6550 Å; in this case, the line is attributed to Hα, and its profile is quite different from SN 2023aew, with a periodic change between a redshift and a blueshift compatible with the periodicity in the light curve of SN 2022jli (Chen et al. 2024). A strong [N II] is more common in SNe IIb from low-mass progenitor stars (Jerkstrand et al. 2015; Fang et al. 2019), and its absence in SNe Ic is a sign that the progenitor star has lost its helium layer. Therefore, the presence of some helium in our spectra is qualitatively compatible with the [N II] identification.

Diverse mechanisms have been suggested for these SNe. Das et al. (2024) consider their sample to be the result of late-time eruptive mass loss: the first, short-lived peak in these objects was argued to be the shock breakout from a dense confined CSM or an extended envelope, followed by the second peak powered by ⁵⁶Ni decay. Similarly, a shock breakout is likely responsible for the first peak of SN 2020bvc (Ho et al. 2020; Izzo et al. 2020). In SN 2022xxf, on the other hand, CSI involving a H- and He-free CSM, seen through late-time narrow peaks on top of broad lines, powers the second peak according to Kuncarayakti et al. (2023). In SN 2019cad, Gutiérrez et al. (2021) argue that the double-peaked structure is the result of a two-component distribution of ⁵⁶Ni (where the outer component could be created by jets), with possibly an additional contribution from a magnetar central engine. The possible magnetar would dominate the second peak. Gomez et al. (2021), conversely, suggest that the first peak of the hybrid SLSN/Ic event SN 2019stc is powered by a magnetar and ⁵⁶Ni and the second by CSI. Finally, Chen et al. (2024) suggest the second, slowly evolving peak of SN 2022jl is powered by the accretion of material from a companion star onto the newly formed compact object. The compact object would be orbiting the companion star on an eccentric orbit that it was kicked into by the explosion, explaining the periodicity of the light curve and the Hα profile.

While the presence of two peaks in it cannot be confirmed, some similarities to the second peak of SN 2023aew can be seen in the SLSN I SN 2019szu (Aamer et al. 2024). An increasing temperature following a plateau is seen (albeit much higher than in SN 2023aew and not conforming to a single blackbody), and nebular lines appear in its spectra very early after the main peak. The SN is, however, much more luminous than SN 2023aew, and its inferred ejecta mass (∼30 M_⊙) much higher than what we obtain from the MOSFiT models (Sect. 4.2). This SN also shows no hint of hydrogen or helium lines. Aamer et al. (2024) explain the plateau as shell-shell collisions from pulsational pair instability (PPI; Woosley 2017), while the main peak is inferred to be dominated by a magnetar.

As none of the diverse objects examined here provides a good analog to SN 2023aew, it is not clear which, if any, of these mechanisms are involved; however, the long separation between the peaks clearly argues against a shock breakout followed by a ⁵⁶Ni peak and against a two-component ⁵⁶Ni distribution, as does the fact that a ⁵⁶Ni model fails to account for the shape of the second peak (see Sect. 4.2). We now discuss other possible physical scenarios for SN 2023aew.

5.2. Physical scenarios

5.2.1. Two unrelated SNe

At first glance, the simplest explanation for SN 2023aew is that each peak is a separate SN: first, a SE-SN spectroscopically similar to a SN IIb, followed by another peculiar SN of Type Ic from a different progenitor system along roughly the same line of sight. There is a $0\stackrel{″}{.}090\pm 0\stackrel{″}{.}053$ (∼1.7σ) difference between the centroids of the SN during the first and second peak (see Appendix C for details), which, taken at face value, would tentatively support this scenario. This discrepancy is not statistically significant, but this is a conceivable, though very unlikely, scenario. The two SNe would have to be separated by only ∼50 pc (projected) in space and four months in time, not to mention their location at the outskirts of the host galaxy, where the star-formation rate and numbers of stars are low (Sect. 3.1).

The Hα line flux caught on the $1\stackrel{″}{.}0$ (∼500-pc) slit in our GTC spectra is on the order of 2 × 10⁻¹⁷ erg s⁻¹ cm⁻². From this, even assuming this Hα emission is concentrated within the 50-pc-radius region where the two SNe would be located, we get an Hα luminosity of L_Hα ∼ 10³⁸ erg s⁻¹ in this region, which corresponds to a star formation rate on the order of 10⁻³ M_⊙ yr⁻¹ (Kennicutt 1998), and is low enough to be powered by only several O stars (Crowther 2013). This, in turn, would translate into a CCSN rate on the order of 10⁻⁵ yr⁻¹ assuming the ratio between the CCSN rate and star formation rate in the nearby universe estimated empirically by Botticella et al. (2012), $\sim 0.01{\text{M}}_{\odot }^{-1}$. In a four-month period, this means a chance probability on the order of only 10⁻¹⁰ for two CCSNe at this location, assuming a Poissonian process. This assumes a constant star formation (and CCSN) rate, however, which can result in an underestimate if the star formation, in reality, occurs in distinct episodes – but at the same time, the star formation may occur more evenly over the 500-pc region covered by the slit, offsetting this effect.

This scenario presents some additional problems besides the low odds of the line-of-sight and timing coincidence: it would not only be a coincidence of two SNe, but of two peculiar SNe. The SN Ic needs to have a low velocity, a slow spectral evolution with a post-peak appearance from the beginning, and a somewhat bumpy light curve that does not fit the ⁵⁶Ni model (Sect. 4.2). Even the SN IIb exhibits a flattening light curve inconsistent with ⁵⁶Co decay, and SN 2003bg, which matches the properties of SN 2023aew at this epoch, shows an abnormally large kinetic energy and ejecta mass (Hamuy et al. 2009; Mazzali et al. 2009). For these reasons, we consider it even more unlikely that SN 2023aew consists of two unrelated SNe.

5.2.2. Two SNe from the same system

A similar scenario that does not require an extremely coincidental chance alignment of unrelated SNe is that of two related SNe. One could conceive a scenario where the two stars in a binary system have very precisely the same mass, resulting in almost simultaneous explosions. Such a scenario could involve either two stars too distant to interact or a common-envelope evolution with a double core (e.g., Fig. 6 of Vigna-Gómez et al. 2018). While roughly equal mass ratios (q) can, in principle, exist in binaries, they are somewhat rare (even q > 0.7 was rarely seen among massive stars by Moe & Di Stefano 2017), and the difference in masses would have to be extremely fine-tuned to produce SNe with a delay of months, as opposed to millennia or more.

If we imagine an equal-mass binary system where the two stars never interact, we would expect the two SNe to be the same spectral type, which seems not to be the case. Although the spectra associated with the second peak seem to retain some relatively weak helium features, they most closely resemble SNe Ic (Sect. 3.3) while the initial spectrum shows strong Hα absorption. It is conceivable that the two progenitor stars never interact with each other, but that one or both has its own closer companion in a hierarchical triple or quadruple star system, resulting in different spectral types despite the same lifetime, but this again requires extreme fine-tuning.

Since most SE-SNe are the result of binary interactions (e.g., Smith et al. 2011), we must consider the effect of Roche-lobe overflow on our hypothetical system. In a typical case where q < 1 and one star evolves before its companion, the mass transfer would increase the mass of the envelope of the companion, leading to rejuvenation and delayed death. In this case, we can obtain one stripped-envelope SN, but not two in rapid succession, and the first event would be expected to be more stripped, not less. On the other hand, in an equal-mass binary where both stars evolve simultaneously, one could recover two He stars. This scenario is invoked by Vigna-Gómez et al. (2018) as a step in one of their channels explaining the observed population of neutron star binaries. The challenge is then to explain the Hα line at −101 d. Although Vigna-Gómez et al. (2018) do not expect to retain hydrogen in the envelope, this could be the result of approximations made in rapid population synthesis codes. As highlighted by Laplace et al. (2020), in detailed models that take into account stellar structure, some hydrogen can remain after a common envelope, which could explain the presence of hydrogen features. The hydrogen could then be blown away by the first SN. It is important to note that these considerations are qualitative and that detailed modeling is required to assess the feasibility of this scenario, which is beyond the scope of this paper. Even if the two-SN scenario is, in principle, feasible, the problem of simultaneous explosions and mass fine-tuning remains.

The problem of fine-tuning could be alleviated if one SN can trigger another in the same system. However, we are not aware of models that can achieve this. Hirai et al. (2018) studied the interaction between ejecta and a close companion star, which only resulted in the bloating and temporary brightening of the companion. In this study, the companions are hydrogen-rich, likely having gained mass from the primary on top of their own remaining hydrogen. One could, in principle, imagine another scenario where the ejecta of the first SN interacts with a stripped companion star instead (possibly after a common envelope phase), but fine-tuning of the mass ratio would still be required. Additionally, since only a small fraction of the kinetic energy of the first SN – itself a small fraction of the binding energy of the core – can affect the secondary star, it is very unlikely that the nuclear burning could be appreciably sped up.

As a side note, whether the two progenitor stars interact or not, and whether the two SNe are separate or one is triggered by the other, we also note that the ejecta of the second SN should interact with that of the first. While the velocities are higher in the −101 d spectrum than later, and thus much of the ejecta from the first SN would be moving too fast to be caught up to by the second SN, some of it would be moving slower due to homologous expansion and could be encountered. Some could possibly even still be located around the progenitor of the second SN, and interaction with this material could conceivably slow down the ejecta of the second SN enough to explain the observed low velocities. While, to our knowledge, there are no published models for this, such an interaction between two massive matter “shells” could produce an extremely strong interaction, likely resulting in a SLSN and an interaction photosphere that would block the line of sight to the inner ejecta of the second SN.

5.2.3. Eruption followed by SN

Numerous SNe have been observed with outbursts and eruptions immediately preceding a SN. Prominent examples of this behavior are the Type IIn SN 2009ip (e.g., Mauerhan et al. 2013; Pastorello et al. 2013) and other similar events (e.g., Pastorello et al. 2018; Fransson et al. 2022), where the precursor eruption typically precedes the main peak by roughly a month. Although some debate remains over whether such transients are true SNe or merely caused by collisions of ejected CSM shells from separate eruptions (Fraser et al. 2013; Pessi et al. 2023), the progenitor stars of SN 2009ip and the similar SN 2016jbu have declined below their pre-SN levels (Brennan et al. 2022; Smith et al. 2022), suggesting a terminal explosion. The observed broad, luminous, and bumpy light curve (Sect. 3.2) could be powered in part by interaction with this ejected material. The final eruption could also have removed (most of) the remaining hydrogen and helium of the progenitor star before its explosion, resulting in the SN Ic spectrum of the second event.

However, 2009ip-like SNe mostly exhibit precursors around r-band absolute magnitudes of −15 mag (Pastorello et al. 2018) with at least one exception, SN 2019zrk, reaching −16.5 mag (Fransson et al. 2022), whereas the peak absolute magnitude of any possible eruption in SN 2023aew must reach a peak of ∼ − 17.9 mag in the i band, which is even somewhat bright for a real SN IIb (Richardson et al. 2014) although similar to that of SN 2003bg (Hamuy et al. 2009). The delay between the eruption and the main peak is also longer by a factor of a few.

SN 2019szu (Aamer et al. 2024) does show an early plateau attributed to a pre-explosion shell-shell collision, reaching an absolute magnitude of −18.7 mag, and this plateau may have lasted a long time and/or been preceded by a brighter peak, similarly to SN 2023aew. In this case, however, the CSM velocity is estimated as 1500 km s⁻¹. We observe absorption minimum (i.e., bulk) velocities of > 10 000 km s⁻¹ in the Hα and Ca II lines in the first spectrum at −101 d. Such velocities are typically only seen in true SNe – Smith et al. (2018) do find velocities in this ballpark in light echoes of the Great Eruption of η Carinae, but most of the mass ejected in this eruption moved at a much lower velocity of ∼600 km s⁻¹. A similar high-velocity component is also seen in SN 2009ip (Pastorello et al. 2013), but again, it is only formed by a small fraction of the total mass. Any putative eruption has to replicate the line profiles of a true SN without the narrower core that should contain the majority of the equivalent width of the lines similar to SN 2009ip and its Great Eruption. Even considering the low resolution of SEDM, the core should be unresolved, on the order of 3000 km s⁻¹ wide.

The first spectrum at −101 d resembles that of a post-peak SN IIb, consistently with the light curve. This and the strength of the absorption lines suggest an ejecta mass and velocity, and hence kinetic energy, indicative of a true SN. This is also compatible with the MOSFiT fit to the first event (Sect. 4.2): the model parameters found from the fit are close to those of SN 2003bg, and the kinetic energy ($\sim \frac{3}{10}{M}_{\mathrm{ej}}{v}_{\mathrm{ej}}^2$) with the median parameter values is ∼2.9 × 10⁵¹ erg. Such an energy is theoretically achievable in an eruption with PPI (Woosley 2017). However, PPI eruptions energetic enough to resemble CCSNe are expected to occur in very massive stars with helium cores of ≳50 M_⊙, with low shell velocities (a few thousand km s⁻¹) more similar to SN 2019szu (Aamer et al. 2024) or SNe IIn, Ibn or Icn.

We also note that the spectra of the second peak seem to exhibit multiple nebular lines starting around +70 d, at the same time as the light curve reaches the luminosity of the plateau that precedes the rebrightening, where lines from that component might re-emerge. These lines include [O I], Mg I], [Ca II]+[O II], and possibly [N II]. [Ca II]+[O II], especially, seems to emerge very early, even before the second peak. This is reasonable if the first peak is in fact a true SN IIb and the lines actually emerge at a later time relative to explosion. The velocities of the lines during the second peak are low compared to normal SNe Ic, but more similar to what one would expect at late times from a SN IIb (e.g., SN 2011dh showed an Hα velocity of ∼ − 12 000 km s⁻¹ at +20 d and absorption line velocities of 4000–7000 km s⁻¹ at +100 d; Ergon et al. 2014). The first peak is red, but the ⁵⁶Ni model in MOSFiT has no problem reproducing the color evolution with a low extinction. Considering all of this, we disfavor the eruption+SN scenario and consider the first peak a real SN.

5.2.4. SN with delayed CSI

If the first peak is powered by a real SN IIb, and the second peak cannot be a real SN Ic from another star, we are left with scenarios where the second peak is instead caused by the delayed input from some mechanism that can inject extra energy into the SN well past the initial explosion. One such mechanism could be CSI, the preferred scenario for SN 2022xxf (Kuncarayakti et al. 2023). CSI is also seen at late times in SNe IIb such as SN 1993J (e.g. Patat et al. 1995), albeit not to the degree that SN 2023aew would require. The CSM could be produced in late-time eruptions. The ejecta mass we obtain from the MOSFiT fit to the first peak is much higher than in either of the CSI models, but since those models assume the second peak is a SN, this may not be a problem. The CSM could, in principle, be created through PPI, but in such a case the progenitor must be very massive compared to typical SN IIb progenitors (Woosley 2017), and most of the ejecta should fall back to create a black hole. The plateau following the SN IIb and the bumps after the second peak suggest the presence of some interaction; a similar flattening has been seen in some interacting SE-SNe (e.g., Kuncarayakti et al. 2022; Ferrari et al. 2024). However, it is unclear if CSI is the dominant power source in SN 2023aew.

Interaction with a dense shell of CSM at a relatively large distance from the progenitor star, to account for the delayed second peak, presents some problems. The onset of the CSI powering the second peak (∼−25 d) would be ∼110 d after the explosion. Assuming a velocity of ∼12 500 km s⁻¹, observed in the first spectrum, for the outer ejecta, the distance to the main CSM shell would be ≳10¹⁶ cm. One problem in this scenario is the blackbody evolution (Sect. 4.1.1). While the temperature at the peak epochs is likely an overestimate because of the spotty NIR light curve interpolation, the dependence of the luminosity on temperature is much stronger than on the radius. Therefore the peak-epoch radius is unlikely to be much larger than 2 × 10¹⁵ cm, and the radius is observed to stay fairly constant for about 50 d before continuing to decline. The photosphere never reaches distances even close to 10¹⁶ cm, where the dense shell should be located; even with a diluted photosphere, an increase in radius should be observed from the onset of CSI, but during the rebrightening phase, only the temperature increases. Instead, this behavior implies an internal mechanism heating the ejecta.

We also see a decrease in absorption velocities from > 10 000 km s⁻¹ at −101 d to ∼6500 km s⁻¹ at −19 d, after which the velocity slowly declines to ∼5500 km s⁻¹ over roughly 160 days. This implies that the matter above the photosphere is located further into the ejecta than at −101 d, while the interaction should be triggered by the ejecta further out and obscure some of the slower material from view. The interaction with a distant shell would also result in boxy emission line profiles (e.g., Patat et al. 1995) that are not observed at this time. Since the SN in this scenario is of Type IIb and contains hydrogen and helium in its spectrum, the obvious implication is that both the outer ejecta of the progenitor star and its CSM must also primarily contain hydrogen and helium; therefore, the boxy emission should be seen in Hα, specifically. No narrow emission lines similar to those in SN 2022xxf (Kuncarayakti et al. 2023) are seen at late times, either, except a faint Hα attributable to the host galaxy, accompanied by [O III] λλ4959, 5007 of similar strength.

CSI can provide an internal heating mechanism if, say, a disk-shaped CSM is engulfed by the ejecta. A disk-interaction model has been proposed for various transients Smith et al. (2015), Andrews & Smith (2018), Pursiainen et al. (2022, 2023). A hypothetical disk located at ∼10¹⁶ cm to account for the delay would have the same problems as a shell at this distance. However, it is possible that the interaction begins somewhat earlier, but most emission from it takes until ∼110 d from the explosion to break out. If the outer ejecta encounters the inner edge of the disk ∼25 d after the explosion, this corresponds to ∼2 × 10¹⁵ cm, the blackbody radius we see at the end of the plateau. The main breakout of emission would, in this scenario, not happen until the outer ejecta becomes optically thin enough. Before this, we would only see the plateau, once the SN itself stops dominating the total emission. The breakout could cause the photosphere to stall as it recedes close to the interaction site, and would cause the increase in observed temperature. The reprocessing site of the CSI energy would move inward in the ejecta over time, creating the absorption lines we see. Nagao et al. (2020) model the light curves of SNe II with disk CSI and find that the observed temperature and luminosity rise very rapidly once the interaction site is uncovered, while being affected only slightly before this happens – but we should in this case also see lines from the interaction region, unless only the ejecta around this region is still optically thick. A possible Hα line (albeit broad and without Hβ) is seen in the spectra, and at late times, it looks somewhat boxy or double-peaked, as might be expected from CSI – but it only appears at +40 d or so, and it is also quite possible that the feature is dominated by [N II] instead. Furthermore, the ejecta of a SN IIb is less massive (∼3 M_⊙ according to MOSFiT) and less opaque (fewer free electrons) than in a SN II, making it more difficult to hide the interaction until ∼110 d.

5.2.5. SN with delayed magnetar engine

Another possible internal heating mechanism for replicating the blackbody evolution – a shrinking or stagnant photosphere but a rising temperature during the rebrightening – could conceivably be a central engine, such as the spin-down of a magnetar (Kasen & Bildsten 2010). The observed evolution might be the result of a combination of a delayed onset of the engine and the shrinking of the photosphere with the decreasing optical depth, which is also needed in the CSI scenario.

The MOSFiT-implemented magnetar model (Nicholl et al. 2017) is unable to fit the temperature evolution of the second peak well (see Sect. 4.2), even though it is the best of the models used here. This manifests as a different rise time and peak epoch in each optical band because of cooling, while in SN 2023aew, the peak is observed to be simultaneous in all bands because of a rising blackbody temperature. The slsn model in MOSFiT, however, assumes the ejecta structure of an infant SN, without 100 d of expansion. The photosphere at the epoch of rebrightening was located at ∼2 × 10¹⁵ cm and shrinking, possibly allowing hotter inner layers to emerge over time and changing the typical temperature evolution of the magnetar-powered SN model. After the second peak, the magnetar model is better equipped to replicate the light curve, although at very late times the light-curve features would somewhat deviate from it (due to, for example, a decreasing trapping fraction). In this scenario, relatively weak CSI could produce the bumps in the light curve.

More importantly, however, to our knowledge, no magnetar-powered model includes a delay between the SN and the onset of the central engine – caused by the spindown of the newborn magnetar through magnetic dipole radiation – let alone a delay of over 100 d (or possibly ∼50, similarly to the disk CSI scenario). Magnetar scenarios, even those that have been invoked for double-peaked SNe such as SN 2005bf (Maeda et al. 2007) and SN 2019cad (Gutiérrez et al. 2021), assume an immediate onset of the magnetic braking. It is not clear if a long delay can be produced in reality, and in the case of SN 2023aew, it is needed. We thus disfavor this scenario.

5.2.6. SN with delayed accretion

A third way to heat the ejecta from within could be fallback accretion onto a newly formed black hole. Such a scenario was considered for SN 2005bf by Maeda et al. (2007), but ultimately disfavored because of the long delay required for the onset of the accretion. Moriya et al. (2018b) also suggest fallback accretion to power the highly energetic SN II OGLE-2014-SN-073: in this case, the accretion would begin immediately, but the peak of the light curve would not be reached until ∼100 d after the explosion because of a large ejecta mass.

In the case of SN 2023aew, the delay before the onset of the main accretion episode (possibly preceded by weaker accretion that causes the flattening light curve before the rebrightening) must be even longer. The velocity of the slow inner ejecta should be fine-tuned to be just below the escape velocity in such a way that falling back onto the black hole and/or creating the accretion disk would take months. A part of the inner ejecta would receive additional energy deposited by the accretion disk wind and be re-accelerated. As a result, the density of the inner ejecta would stay higher than expected from ∼100 d of homologous expansion, and it would be possible to have a combination of nebular lines from the outer ejecta and a photospheric spectrum from the inner ejecta. We also note that the line profiles of forbidden oxygen and calcium lines are usually not both double-peaked because of their different distribution in the ejecta (Fang et al. 2024), but in SN 2023aew neither feature is single-peaked at late times (Fig. 8). This may suggest that the forbidden lines are only emitted by the part of the ejecta that did not fall back onto the compact object, perhaps supporting this scenario. If the accretion does not proceed smoothly, bumps in the light curve (and the fast decline compared to the MOSFiTfallback model; see Sect. 4.2) could be explained as well, but an additional contribution from CSI is possible. We encourage detailed modeling of this scenario to determine its plausibility.

The total radiated energy during the second peak, within the observed timescale and wavelength range, is ∼4 × 10⁴⁹ erg, and the peak luminosity is ∼8 × 10⁴² erg s⁻¹ (Sect. 4.1.2). The outflow energy that can be extracted from accretion in a disk wind is E_acc = ϵMc² (and the corresponding energy injection rate is $\frac{\mathrm{d}{E}_{\mathrm{acc}}}{\mathrm{d}t}=ϵ\frac{\mathrm{d}M}{\mathrm{d}t}{c}^2$), where ϵ ∼ 10⁻³ is the radiative efficiency of the accretion (Dexter & Kasen 2013). This corresponds to the luminosity of a directly observed accretion disk, which in our scenario is hidden by the SN ejecta. Assuming that this energy is wholly deposited into the ejecta and roughly 1% can escape as radiation, we can estimate that in order to power the observed light curve, the accreted mass needs to be on the order of a few M_⊙ (note that this is only the part of the fallback that forms the accretion disk and does not fall directly into the black hole). The accretion rate $\frac{\mathrm{d}M}{\mathrm{d}t}$ at the peak would be a few × 10⁻² M_⊙ d⁻¹ with the same assumptions. As the exact processes are unclear, these numbers are simply order-of-magnitude estimates. In this scenario, the fallback of the inner ejecta would lessen the amount of ⁵⁶Ni available for powering the first event, but since the inner ejecta is probably aspherically distributed (Fang et al. 2024), this may not present a problem. The progenitor star would likely be much more massive than those of normal SNe IIb – it would need to have at least several M_⊙ of ejecta divided between the accreted inner ejecta, the accelerated inner ejecta and the outer ejecta responsible for the first peak.

The Eddington luminosity of an accreting compact object is $L_{\mathrm{Edd}}\approx 1.2\times {10}^{38}\frac{M_{\mathrm{CO}}}{{\mathrm{M}}_{\odot }}$ erg s⁻¹, where M_CO is the mass of the object. For a 10-M_⊙ black hole, one thus obtains L_Edd ≈ 1.2 × 10³⁹ erg s⁻¹, roughly a factor of 6700 below even the observed bolometric luminosity at the peak, let alone the total required outflow energy. A similarly super-Eddington luminosity was noted for the Type Ic SN 2022jli by Chen et al. (2024) when using a neutron star as the compact object; SN 2023aew was somewhat more luminous than SN 2022jli, requiring a larger compact-object mass or an even more dramatically super-Eddington accretion. The mechanisms of such accretion are unclear (Brightman et al. 2019), but examples are known, from slightly super-Eddington supermassive black holes through some X-ray binaries (one of which exceeds its L_Edd by a factor of ∼500; Israel et al. 2017) to accretion-powered long gamma-ray bursts (LGRBs). In LGRBs, accretion rates onto stellar-mass black holes can momentarily range from 0.01 to several M_⊙ s⁻¹ (e.g., Kumar & Zhang 2015), orders of magnitude above our estimate for SN 2023aew. The SNe accompanying such events tend to be of Type Ic-BL (e.g., SN 1998bw; Galama et al. 1998). We see no very broad lines indicating an outflow similar to the SNe accompanying LGRBs (indeed, the velocity during the second peak is lower than in most SNe Ic); instead, in SN 2023aew any Ic-BL-like outflow or jet would either be directed outside our line of sight or be too weak to accelerate the (inner) ejecta to such velocities.

A compact object kicked in the direction of a secondary star, similarly to the onset of the second peak in SN 2022jli (Chen et al. 2024), could also trigger accretion, possibly falling inside the envelope of the secondary and causing an outflow of hydrogen-poor material without the periodicity of SN 2022jli. We note that the bumps in the SN 2023aew light curve could have a ∼75-day periodicity, but with only two bumps this is difficult to say for certain, and such a period would imply a considerably larger orbital separation from the companion than in SN 2022jli, where the period was ∼12 d, resulting in weaker accretion except near the pericenter of an eccentric orbit.

6. Conclusions

We spectroscopically and photometrically followed up the peculiar, double-peaked SE-SN, SN 2023aew, in the optical and NIR. This object first shows behavior similar to a SN IIb – albeit somewhat luminous (M_TESS, peak = −17.88 ± 0.12 mag), red and with weak helium lines – followed by a plateau ∼50 d after the first peak and a rebrightening at ∼100 d after the first peak. The second peak resembles a slowly evolving SN Ic with a bumpy light curve and a peak absolute magnitude (M_r = −18.75 ± 0.04 mag) similar to SNe Ic-BL. We analyzed the data, performed fits to the spectral energy distributions and light curves and compared the behavior of SN 2023aew to other SNe Ic and double-peaked SE-SNe. We then explored various physical scenarios to attempt to explain the anomalous behavior of the object. Our conclusions are as follows:

• SN 2023aew, among the double-peaked SE-SNe examined here, presents a unique combination of properties, with very widely separated peaks and a slow spectroscopic evolution. The first peak before the flattening is in many ways (spectrum, ejecta mass, luminosity) similar to the Type IIb SN 2003bg. Compared to normal SNe Ic, the second peak exhibits a low expansion velocity in addition to the slow spectroscopic and photometric evolution.

• The second peak is difficult to fit with a ⁵⁶Ni decay or magnetar central engine model. A part of this difficulty is the simultaneous peak in all optical bands, caused by an increasing rather than decreasing temperature – which, furthermore, only reaches ∼8000 K at peak – during the rebrightening. Meanwhile, the blackbody radius slowly shrinks during almost the entire observed light curve.

• The simplest explanation for the two widely separated peaks, two SE-SNe coincident in both time (∼120 d) and space (≲50 pc projected), is extremely unlikely, considering not only the low star-formation rate in the vicinity of SN 2023aew, but also its peculiarity. Meanwhile, a scenario of two SE-SNe from the same system avoids this problem but requires extremely precise fine-tuning and has difficulty explaining the properties of the two SNe. Therefore, we do not consider a two-SN scenario plausible.

• The line velocities, luminosity and kinetic energy of the first peak are difficult to explain if it is a pre-SN eruption instead of a true SN, whereas these properties and the time scale are compatible with a SN IIb. We thus consider the first peak to be associated with a SN, followed by delayed energy input by some mechanism.

• This mechanism is unlikely to be interaction with a distant (≳10¹⁶ cm) CSM shell, for several reasons. No spectroscopic signatures of interaction are seen; hiding the interaction within a photosphere at such a distance would be difficult; and absorption lines would be formed in the outer ejecta, which is not the case.

• The blackbody parameter evolution instead indicates a source of heating from within the ejecta. This could be emission from embedded CSI with a disk of CSM engulfed by the outer ejecta. The main emission could stay hidden until the outer ejecta becomes optically thin enough. However, the difficulty of hiding the interaction until the second peak and a lack of CSM lines at this time argue against this scenario.

• If the internal source is a magnetar central engine instead, the onset of the magnetar spindown has to be delayed by months through an as-yet unknown and unprecedented mechanism.

• Fallback accretion onto a newly formed compact object could launch a disk wind that re-accelerates the inner ejecta. The long delay before the onset of the accretion requires a fine-tuned, low inner ejecta velocity.

While we can therefore exclude some of the scenarios we have considered, such as two-SN scenarios and pre-SN eruptions, we cannot definitively establish the power source of the second, brighter peak. CSI with a distant shell and some internal heating mechanisms present problems with the time scale and/or the spectroscopic evolution, but a scenario where the second peak is powered by accretion onto a black hole may be the least problematic. We encourage detailed modeling studies of SN 2023aew, outside the scope of this paper, in order to shed more light on this unprecedented SN.

Acknowledgments

We thank the anonymous referee for their comments on the manuscript, and Luc Dessart for helpful discussions. HK was funded by the Research Council of Finland projects 324504, 328898, and 353019. TN acknowledges support from the Research Council of Finland projects 324504 and 328898. RK acknowledges support via Research Council of Finland (grant 340613). MF is supported by a Royal Society – Science Foundation Ireland University Research Fellowship. HFS is supported by the Eric and Wendy Schmidt A.I. in Science Fellowship. KM acknowledges support from JSPS KAKENHI grant (JP20H04737) and the JSPS Open Partnership Bilateral Joint Research Project between Japan and Finland (JPJSBP120229923). MDS is funded by the Independent Research Fund Denmark (IRFD, grant number 10.46540/2032-00022B. CPG acknowledges financial support from the Secretary of Universities and Research (Government of Catalonia) and by the Horizon 2020 Research and Innovation Programme of the European Union under the Marie Skłodowska-Curie and the Beatriu de Pinós 2021 BP 00168 programme, from the Spanish Ministerio de Ciencia e Innovación (MCIN) and the Agencia Estatal de Investigación (AEI) 10.13039/501100011033 under the PID2020-115253GA-I00 HOSTFLOWS project, and the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. Y-ZC is supported by the National Natural Science Foundation of China (NSFC, Grant No. 12303054), the Yunnan Fundamental Research Projects (Grant No. 202401AU070063) and the International Centre of Supernovae, Yunnan Key Laboratory (No. 202302AN360001). MP acknowledges support from a UK Research and Innovation Fellowship (MR/T020784/1). PC acknowledges support via Research Council of Finland grant 340613. SM acknowledges support from the Research Council of Finland project 350458. NER, AP, SB, AR and GV acknowledge support from the PRIN-INAF 2022, ‘Shedding light on the nature of gap transients: from the observations to the models’. NER also acknowledges partial support from MIUR, PRIN 2017 (grant 20179ZF5KS). AR also acknowledges financial support from the INAF GRAWITA Large Program Grant. We acknowledge ESA Gaia, DPAC and the Photometric Science Alerts Team (http://gsaweb.ast.cam.ac.uk/alerts). The Starlink software (Currie et al. 2014) is currently supported by the East Asian Observatory. Partially based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. The data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOT. This study is based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. Access to the Las Cumbres Observatory was made possible via an allocation by OPTICON (program 22A/012, PI Stritzinger). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004719. The GROWTH India Telescope (GIT) is a 70-cm telescope with a 0.7-degree field of view, set up by the Indian Institute of Astrophysics (IIA) and the Indian Institute of Technology Bombay (IITB) with funding from Indo-US Science and Technology Forum and the Science and Engineering Research Board, Department of Science and Technology, Government of India. It is located at the Indian Astronomical Observatory (IAO, Hanle). We acknowledge funding by the IITB alumni batch of 1994, which partially supports the operation of the telescope. Telescope technical details are available at https://sites.google.com/view/growthindia/. Operation of Nayuta telescope and NIC is partially supported by the Optical and Infrared Synergetic Telescopes for Education and Research (OISTER) program funded by the MEXT of Japan. The spectra taken by the Seimei telescope are obtained under the KASTOR (Kanata And Seimei Transient Observation Regime) campaign. The Seimei telescope at the Okayama Observatory is jointly operated by Kyoto University and the Astronomical Observatory of Japan (NAOJ), with assistance provided by the OISTER program.

References

Appendix A: Supplementary tables

Appendix B: MOSFiT corner plots

Fig. B.1. MOSFiT corner plot for the 56Ni decay fit of the first peak.

Fig. B.2. MOSFiT corner plot for the 56Ni decay fit of the second peak.

Fig. B.3. MOSFiT corner plot for the magnetar spindown fit of the second peak.

Fig. B.4. MOSFiT corner plot for the CSI fit of the second peak when s = 0.

Fig. B.5. MOSFiT corner plot for the CSI fit of the second peak when s = 2.

Fig. B.6. MOSFiT corner plot for the fallback fit of the second peak.

Appendix C: Coordinate discrepancy

Considering the differences between the spectra of the first and second peak and the rapidity of the rise to the second peak, the two peaks of SN 2023aew might, in fact, conceivably be two SNe on slightly different lines of sight, although such a spatial and temporal coincidence would be extremely unlikely, especially taking into account the faintness and low star formation rate of the host galaxy (see Sect. 5). In order to check whether this is the case, we used the coordinates of Gaia23ate in the Gaia Alerts Index (α₁ = 17:40:51.37; δ₁ = 66:12:22.75). These correspond to the first scan of the first peak. The uncertainty of Gaia coordinates is 0$\stackrel{″}{.}$055 per scan (Hodgkin et al. 2021). We also searched for highly precise Gaia coordinates of stars within the field, using 26 stars within 2′ of the SN with astrometric uncertainties between 0.02 and 3 mas as reported in the Gaia Archive.³¹ This uncertainty is negligible compared to that of the position of SN 2023aew. We then used the IRAF tasks ccmap, cctran and ccsetwcs to calibrate the astrometry in our ALFOSC and NOTCam images from 2023 Apr 26, 2023 Jul 08 and 2023 May 30; the root-mean-square (RMS) errors from these fits were between 0$\stackrel{″}{.}$03 and 0$\stackrel{″}{.}$06. After rejecting outliers, between 15 and 20 stars visible in the field were used in each fit.

From our NOT images, we measured the average coordinates of the second peak as $\overline{{\alpha }_2}$ = 17:40:51.356; $\overline{{\delta }_2}$ = 66:12:22.820. The uncertainty of any discrepancy is $\sqrt{{\sigma }_1^2+{\sigma }_2^2}$, where σ₁ is the Gaia uncertainty of the first peak coordinates (0$\stackrel{″}{.}$055) and σ₂ the uncertainty of the average NOT centroid measurement, ${\sigma }_2=\sqrt{\mathrm{\Sigma }{\sigma }_i^2}/N$, where N is the number of measurements (9). Each individual measurement error σ_i is calculated as

$$\begin{array}{c}\hfill {\sigma }_i=d_i(\frac{{\sigma }_{i,\mathrm{RA}}}{\mathrm{\Delta }{\mathrm{RA}}_i}+\frac{{\sigma }_{i,\mathrm{Dec}}}{\mathrm{\Delta }{\mathrm{Dec}}_i}),\end{array}$$(C.1)

where d is the angular distance between the peaks in measurement i, ΔRA is the angular distance in the right ascension direction between the peaks (i.e., not the difference in α directly) and ΔDec the declination difference, σ_RA the RMS error of the astrometry in α and σ_Dec in δ. This is because $d=\sqrt{\mathrm{\Delta }{\mathrm{RA}}^2+\mathrm{\Delta }{\mathrm{Dec}}^2}$. We then obtain an offset of $0\stackrel{″}{.}104\pm 0\stackrel{″}{.}080$, with a significance of 1.3σ.

Additionally, we used the gr-band centroid coordinates reported for the ZTF photometry by ALeRCE, assuming that the astrometry is consistent between each epoch (the accuracy of the absolute astrometry is irrelevant here). We took the averages of the α and δ coordinates of all epochs before MJD = 60043 (17 measurements) to represent the first peak and those of all epochs between MJD = 60044 and MJD = 60156 (76 measurements) to represent the second peak. The uncertainties in α and δ were taken to be the standard errors of the mean. The resulting coordinate discrepancy between the two peaks is $0\stackrel{″}{.}078\pm 0\stackrel{″}{.}068$, also barely over 1σ, but in the same direction (negative α, positive δ) as the Gaia-NOT discrepancy.

Finally, combining the α and δ offsets from the Gaia-NOT comparison and the ZTF measurements, we obtain a total offset of $0\stackrel{″}{.}090\pm 0\stackrel{″}{.}053$ if the ZTF and NOT offsets are given equal weight. The significance of this discrepancy is only ∼1.7σ – as such, even if there is any coordinate offset between the two objects, it is too small to significantly detect in the data available to us. There is also a possibility of underestimated errors because of a systematic error in the Gaia coordinates of the first peak, as the offset from the NOT images relies on this measurement. Such an issue could further decrease the significance.

All Tables

All Figures

	Fig. 1. Three-color (blue: g, green: r, red: i) cutout of the field of SN 2023aew, observed using NOT/ALFOSC on MJD = 60134.0. Nuclei of host galaxy candidates, Galaxy A and Galaxy B, are labeled. These galaxies are located at z = 0.30 behind the true host galaxy, Galaxy C (see Sect. 3.1).
In the text

	Fig. 2. Comparison of the SN 2023aew classification spectrum (red; Wise et al. 2023) and SN 2003bg (black; Hamuy et al. 2009), a SN IIb, at +20 d. The resemblance is fairly strong, but the SN 2023aew spectrum is quite noisy (uncertainties are shown as shaded regions). Telluric absorption lines (⊕) are marked with vertical gray lines.
In the text

Fig. 3. NOT/ALFOSC spectra of Galaxies A (black) and B (blue) and the GTC/OSIRIS spectrum of SN 2023aew at +131 d showing narrow host lines (purple). Savitzky–Golay smoothing has been applied to the original ALFOSC spectra (gray). Identified lines of [OII], [OIII], and H at z = 0.30 are plotted as black dashed vertical lines. The same lines are also plotted at the redshift of SN 2023aew, z = 0.025, as red dotted vertical lines. No lines are detected at this redshift in Galaxies A and B.

In the text

Fig. 4. NUV, optical, and NIR light curve of SN 2023aew. G-band errors are not reported at Gaia Alerts, and only the distance uncertainty is included. The shaded regions correspond to the 68% uncertainties of our GP fits. The decay rate of 56Co is indicated with a dashed black line. Our optical spectral coverage is indicated with gray vertical lines and our NIR spectrum with a black vertical line. The epoch of the classification spectrum (Wise et al. 2023) is marked with a blue vertical line. All magnitudes have been corrected for Galactic extinction. Phases refer to the second r-band peak (MJD = 60074.9).

In the text

	Fig. 5. Measured colors of SN 2023aew in various bands when roughly simultaneous (Δt < 0.1 d) epochs are available (points) and subtractions between our GP fits of individual light curves when simultaneous photometry points at the different bands are scarce (lines and shaded regions).
In the text

Fig. 6. Comparison of the light curve of SN 2023aew (points) to other SE-SNe (lines). The first peak at first resembles SN 2003bg (dashed purple line; Hamuy et al. 2009), until the flattening around −100 d. The second peak has a luminosity similar to the Type Ic-BL SN 1998bw (Galama et al. 1998; McKenzie & Schaefer 1999) and LSNe (Gomez et al. 2022) (solid lines), but a much slower evolution. Of the double-peaked SE-SNe (DP; dotted lines) pictured here – SNe 2019cad (Gutiérrez et al. 2021), 2020bvc (Ho et al. 2020; Izzo et al. 2020), 2022xxf (Kuncarayakti et al. 2023), and 2022jli (Chen et al. 2024) – only SN 2020bvc resembles the late-time light curve of SN 2023aew, and the separation between peaks is much shorter in all of them than in SN 2023aew.

In the text

Fig. 7. Spectral evolution of the first (blue; Wise et al. 2023) and second peak of SN 2023aew. Different instruments are color-coded, and epochs refer to the second peak. Telluric absorption lines (⊕) are marked with solid light gray lines. Dash-dotted lines are at zero velocity (of the stronger line in case of a close doublet); long dashes at −2500 km s−1; short dashes at −6500 km s−1; and dotted lines at −8500 km s−1. When Savitzky–Golay smoothing has been applied, the original spectrum is plotted in gray.

In the text

Fig. 8. Left: evolution of the Na Iλλ5890,5986 line (dashed gray lines), possibly blended with He Iλ5876. The approximate evolution of the absorption minimum from a Gaussian fit after −19 d is marked with a solid gray line. The velocity evolution is slow during the second peak and the subsequent decline. Middle: evolution of the putative [O I] λλ6300,6364 (dashed red lines) and [N II] λλ6548,6584 (dashed gray lines) doublets after the second peak. Both start to appear between +41 and +71 d and strengthen over time. Right: evolution of the [Ca II] λλ7291,7323 (dashed gray lines) + [O II] λλ7320,7330 (dashed red lines) feature. Zero velocity refers to the stronger line of each doublet. All forbidden features shown here are double-peaked at late times.

In the text

	Fig. 9. TNG/NICS NIR spectrum of SN 2023aew at −10 d. Possible features of H, He, Mg and C are seen, but the spectrum is noisy and atypical for SE-SNe. Line velocities are as in Fig. 7, with the addition of −12 500 km s−1 as a solid line. Savitzky–Golay smoothing has been applied.
In the text

	Fig. 10. Spectra of SN 2023aew compared to well-observed ordinary SNe Ic. While neither SN 1994I (thick blue line) nor SN 2007gr (dashed red line) matches the early-time features of SN 2023aew, the latter does resemble SN 2007gr at later epochs. The feature at ∼6500 Å is not seen in late-time spectra of either SN 1994I or SN 2007gr. Overall the evolution of SN 2023aew is relatively slow.
In the text

	Fig. 11. Evolution of the Na Iλλ5890,5896 absorption minimum velocity in SN 2023aew (black) after the rebrightening, compared to SNe 1994I (blue) and 2007gr (red). Uncertainties are estimated at 5%. For SN 2023aew, the phase refers to the second peak.
In the text

Fig. 12. Spectra of SN 2023aew (black) compared to double-peaked SE-SNe SN 2019cad (purple), SN 2022jli (thick green line) and SN 2022xxf (dashed brown line). All epochs refer to the second peak. SN 2023aew evolves slower than SNe 2019cad or 2022xxf and exhibits lower velocities, but an excellent match is seen with the latest spectrum of SN 2019cad. The 6500 Å feature, in particular, is replicated. The match with SN 2022xxf similarly improves at late times, although this feature does not appear. SN 2022jli provides a good match at ∼30–40 days and evolves slower than other comparison objects, but lacks a strong [O I] doublet at late times.

In the text

	Fig. 13. Blackbody fits (lines) to the observed SEDs (points) at roughly 10-day intervals. NUV points (light gray) at +2 and +10 d clearly deviate from a blackbody because of line blanketing, and were ignored in the fits.
In the text

	Fig. 14. Evolution of the blackbody parameters over time. The temperature from the fit rises during the rebrightening from a preceding plateau and is highest at the peak – albeit likely slightly overestimated – then declines to ∼5500 K and stays constant until ≳120 d. The radius slowly declines almost throughout the observed evolution.
In the text

Fig. 15. SuperBol (Nicholl 2018) pseudo-bolometric light curve of the second peak integrated over observed bands (black open circles), and a bolometric light curve based on blackbody fits (red filled circles), highly uncertain before the rebrightening. The difference is notable only around maximum light. The late-time decline is inconsistent with 56Co decay (solid black line). Pseudo-bolometric light curves of SNe 1998bw (Patat et al. 2001), 2019cad (Gutiérrez et al. 2021), 2022jli (Moore et al. 2023), and 2022xxf (Kuncarayakti et al. 2023) are included for comparison.

In the text

Fig. 16. MOSFiT light curve fits to the first and second peak of SN 2023aew using the different models described in the text. The flat part of the light curve after MJD = 60000 was ignored in the first-peak fit, as the 56Ni model does not include such features. All models have trouble reproducing the simultaneity of the second peak in all optical bands, specifically failing to reproduce the u-band peak, but by eye and by MOSFiT score, the magnetar model provides the best fit.

In the text

	Fig. B.1. MOSFiT corner plot for the 56Ni decay fit of the first peak.
In the text

	Fig. B.2. MOSFiT corner plot for the 56Ni decay fit of the second peak.
In the text

	Fig. B.3. MOSFiT corner plot for the magnetar spindown fit of the second peak.
In the text

	Fig. B.4. MOSFiT corner plot for the CSI fit of the second peak when s = 0.
In the text

	Fig. B.5. MOSFiT corner plot for the CSI fit of the second peak when s = 2.
In the text

	Fig. B.6. MOSFiT corner plot for the fallback fit of the second peak.
In the text