© The Authors 2026

1. Introduction

Superluminous supernovae (SLSNe) are a recently recognized class of astronomical transients whose optical luminosities exceed those of typical core-collapse supernovae (SNe) by a factor of 10–100. Reaching peak magnitude typically 40 days after the explosion, some SLSNe have reached a u band absolute magnitude brighter than −22 mag, which translates to luminosities in the range of 10⁴⁴–10⁴⁵ erg s⁻¹. Quimby et al. (2011) noted that the late-time decay rates of these particular SNe are inconsistent with a radioactive decay scenario and require a late deposition of a large amount of energy to explain their light curves. As for the early classification of SNe, SLSNe are divided into two main categories based on their spectral properties; namely, the absence or presence of hydrogen lines, corresponding, respectively, to type I (hydrogen-poor) and type II SLSNe (hydrogen-rich). Type I SLSNe are the most studied class and there are now over 260 SLSNe reported in the catalog of Gomez et al. (2024).

What makes SLSNe so different from regular core-collapse SNe is still being debated. There are mainly four different energy sources being considered to explain the high peak luminosity of SLSNe: ejecta fallback accretion onto a black hole, radioactive decay of ⁵⁶Ni, circumstellar interaction, and a magnetar wind nebula (see, e.g., Moriya et al. 2018b; Gal-Yam 2019; Inserra 2019; Nicholl 2021; Moriya 2026, for a review). In the last of these scenarios, a neutron star with extreme properties (millisecond period and ∼10¹⁴–10¹⁵ G magnetic field; a magnetar) will convert its rotational energy into a wind of electrons and positrons. This magnetar wind nebula will produce high-energy nonthermal radiation (via inverse Compton (IC) and synchrotron) that upon interaction with the surrounding ejecta will be reprocessed into lower-energy photons (optical/UV) and power the SLSNe. In theoretical models aiming to predict the optical light curves from engine-powered SLSN (e.g., Dessart et al. 2012; Chatzopoulos et al. 2013; Metzger et al. 2015), it is generally assumed that the thermalization efficiency is 100% (all high-energy radiation is reprocessed by the ejecta). While this is certainly a viable assumption in the early phase where the ejecta are opaque, this efficiency will depart from 100% as the ejecta expands and some γ-rays start leaking. Therefore models predicting the light curve of escaping γ-rays need to treat this efficiency carefully in a time-dependent manner, requiring the computation of the diffusion and escape of the photons through the ejecta (for more details, see Vurm & Metzger 2021). Note that while the aforementioned model is focused on a magnetar as central engine, most of the general principles outlined in the model remain valid for other central engines capable of injecting a pool of nonthermal radiation. An alternative scenario is the circumstellar medium (CSM) interaction, which could also account for most thermal observables (e.g., see Chatzopoulos et al. 2013). It proposes that the extreme luminosity arises from the conversion of the SN ejecta’s kinetic energy into radiation. When the fast-moving ejecta collides with a dense shell of CSM, a strong shock is created that transforms kinetic energy into thermal energy radiating mostly in the X-ray band. This high-energy emission is absorbed by the optically thick CSM and then reemitted as optical/UV light. The slow diffusion of the photons through the CSM broadens the event’s duration, accounting for the observed long-lasting light curves. Observations of some type I SLSNe such as SN 2017egm have revealed irregularities in their light curves, including post-peak bumps, which require more complex multiple CSM shell scenarios (Lin et al. 2023). Other scenarios such as radioactive decay or accretion could play a role in some objects but are for now not thought to be the dominant mechanism. For example, based on their SLSNe catalog, Gomez et al. (2024) concluded that only a small fraction of the sample can accommodate a significant radioactive decay component. For the accretion scenario, Moriya et al. (2018a) studied the required accreted mass to power the SLSNe light curve and concluded that in most cases the mass is greater than the ejecta mass.

The search for SLSNe in γ-rays is more recent and only upper limits were reported up to the claimed γ-ray detection of SN 2017egm by Li et al. (2026, see our Sects. 1 and 3.4). Using Fermi-LAT data, Renault-Tinacci et al. (2018) presented individual and stacked analysis of a sample of 45 SLSNe and only upper limits were reported. However their samples ranged from 2008 to 2015 and the nearest object tested (CSS140222) was at a distance of 145 Mpc with relatively poor multiwavelength coverage. All other sources investigated lie at a redshift of z > 0.1 (d > 400 Mpc). Naturally SN 2017egm was not present in the sample.

Regarding Type IIn (i.e., with narrow emission hydrogen lines, conveying strong CSM interaction), individual studies (e.g., SN 2009ip, Margutti et al. 2014) or systematic searches have been carried out, as in Ackermann et al. (2015) on 147 SNe Type IIn evolving in dense CSM in a one year time window. No significant excess was detected, including the most promising source SN 2010jl (at 49 Mpc; see also Murase et al. 2019). Another γ-ray study was carried out by Prokhorov et al. (2021), who did a systematic search with Fermi-LAT from a large sample of more than 55 000 SNe candidates from the open SN catalog (not only SLSNe). Their analysis used a sliding time window technique in an aperture photometry mode that allowed for the exploration of a large number of candidates but that is less sensitive than a full forward folding likelihood analysis such as the one that we present in this work for a small sample of SLSNe. A few candidate detections were reported, including variable emission spatially coincident with the SN iPTF 14hls. However, a detailed study by Yuan et al. (2018) reported that the γ-ray emission lies within 0.045° of a blazar, which offers an alternative scenario for the high-energy signal.

Acharyya et al. (2023) present a γ-ray follow-up with Fermi-LAT and VERITAS telescopes of the nearby SLSNe SN 2015bn and SN 2017egm, and only upper limits are reported. We note that in a 6 month time bin for SN 2017egm with Fermi-LAT data an interesting hint of a signal (test statistic TS of 10.1) is reported but not considered further.

SN 2017egm was first reported by the Gaia telescope as Gaia17biu on May 23, 2017. At a redshift of z = 0.03 (135 Mpc¹), it was the closest type I SLSNe at the time of detection². In the spectral analysis, the source revealed no H emission (Type I SLSN) and He I emission lines classifying the source as member of the small but growing class of helium-rich SLSN-Ib (Zhu et al. 2023). The host galaxy of SN 2017egm (NGC 3191) is a massive, metal-rich spiral galaxy, in contrast with most Type I SLSNe, which are mostly found in metal-poor dwarf galaxies (Nicholl et al. 2017). While the early-time light curve of the first months can be reproduced with only a magnetar model in Nicholl et al. (2017), the late-time observations revealed multiple bumps starting 200 days after the SN that cannot be explained by either a magnetar model alone or a one-zone CSM interaction (Zhu et al. 2023; Lin et al. 2023). Although temporal variation in the properties of the magnetar and in the opacity could be invoked to modulate the observed flux in a magnetar-dominated model (see Hosseinzadeh et al. 2022), the scenario favored by Lin et al. (2023) to explain the bumpy features is the multiple CSM shell interaction plus radioactive decay model. This multiple-shell scenario could result from pulsational mass ejections from the SN progenitor due to the electron-positron pair instability mechanism that ejects large quantities of matter without totally disrupting the progenitor (Woosley et al. 2007; Chatzopoulos & Wheeler 2012; Woosley 2017; Renzo et al. 2020).

At high radio frequencies (5–34 GHz), non-detections were reported in the 34–74 day time interval after the explosion (Coppejans et al. 2018). In the X-ray band, Chandra follow-up observations obtained throughout the first year post-explosion were reported in Zhu et al. (2023) and yielded the most stringent upper limits of any SLSN-I to date. At higher energies, the recently reported detection of giga-electronvolt γ-ray emission from SN 2017egm by Li et al. (2026) with Fermi-LAT adds another ingredient to the understanding of this complex system.

This triggered our interest to revisit a sample of nearby (< 200 Mpc) sources based on the curated lists of Type I SLSNe by Gomez et al. (2024) and type II SLSNe by Pessi et al. (2025). In Sect. 2 we define a list of interesting SLSNe candidates for γ-ray follow-up with Fermi-LAT and present the light curves for the 16 yrs of data available with Fermi-LAT. In Sect. 3 we carry out a detailed analysis of the most interesting target SN 2017egm and compare our results with those of Li et al. (2026), while we confront the temporal and spectral properties with the models of Vurm & Metzger (2021) in Sect. 5. Predictions for the detection of such a target at tera-electronvolt γ-rays with the Cherenkov telescope CTAO are then explored in Sect. 7 and alternative multiwavelength counterparts to the γ-ray signal are discussed in Sect. 4.

2. Sample selection and data analysis

2.1. SLSN sample selection

The classification of a transient source as a SLSN is non trivial and can evolve over months to year timescales after the SN discovery as more data become available. In particular, a critical parameter to claim the SN to be luminous or super-luminous is the distance to the source. Depending on the quality of the spectral data obtained, the association with a host galaxy and the redshift estimation can evolve with time but not always end up being updated in online databases such as WISeREP³ or TNS⁴. In addition, some sources can be later reclassified in another category such as tidal disruption events or other classes as discussed in Appendix D of Gomez et al. (2024) and Sect. 6 about contaminants in Pessi et al. (2025). Therefore, we decided to base our γ-ray follow-up on curated lists of events to ensure a robust distance estimation and classification. We selected the Type I SLSNe (H-poor) catalog of Gomez et al. (2024), which provides a verified sample of 262 SLSNe reported through the end of 2022. For type II SLSNe (H-rich), we used the catalog of Pessi et al. (2025), which comprise a sample of 107 objects. Combining both catalogs with a distance threshold of 200 Mpc, we end up with the list of 6 targets presented in Table 1. Such a threshold is motivated by the Fermi-LAT sensitivity: a SLSN at 200 Mpc with a γ-ray luminosity similar to that in the optical band (∼10⁴⁴ erg s⁻¹) would produce an energy flux comparable to the LAT 2σ sensitivity level for 1 month of exposure (obtained by scaling the upper limit derived in the SN 2023ixf Fermi-LAT study; Martí-Devesa et al. 2024).

2.2. Optical properties of the sample

To compare the optical properties of our sample in a self-consistent way, we retrieved the publicly available light curves of each object from the literature and public archives to reconstruct the r band absolute magnitudes and the pseudo-bolometric luminosities for each SLSN with a set of common bands. For this, we used data from various sources: Zhu et al. (2023) for SN 2017egm; Chen et al. (2021) for SN 2018bsz; Gomez et al. (2024) for SN 2019ieh; Tinyanont et al. (2023) for SN 2020wtn; Brennan et al. (2024) and Pessi et al. (2025) for SN 2021adxl; and Pessi et al. (2025) for SN 2022mma. Data from the ZTF data release were also used⁵. SN 2019ieh is technically a luminous SN I (Gomez et al. 2022) but we include it here as it is part of the sample of SLSN I of Gomez et al. (2024).

We note that the only photometric bands that are commonly available for all objects of our sample are g, r and i. Thus, we constructed pseudo-bolometric light curves considering only these bands following the second method presented by Pessi et al. (2025), which considers the integration of the spectral energy distribution (SED) over the interpolated gri band light curves with respect to r band peak, considering extrapolations to both the UV and NIR by fitting a blackbody to consider the missing flux. In the UV flux, we extrapolate the blackbody fit from 0 Å to our g band and in the NIR flux, i band to 25 000 Å. To compare the light curves, we aligned them to their r band peak epoch as reported in Chen et al. (2021) for SN 2018bsz (MJD 58267.5), Gomez et al. (2024) for SN 2019ieh (MJD 58671.6), Tinyanont et al. (2023) for SN 2020wnt (59213.75 MJD), Brennan et al. (2024) for SN 2021adxl⁶ (MJD 59521), Pessi et al. (2025) for SN 2022mma (MJD 59790.76). For SN 2017egm the reported peak has been calculated on the g band by Bose et al. (2018). Using Gaussian process interpolation we obtained the r band peak epoch of SN 2017egm to be MJD 57926.4 ± 1.

Our pseudo-bolometric light curves differ from those available in the literature for each object, either because different cosmology parameters are considered or additional bands are included. Although the pseudo-bolometric luminosities obtained in this work only provide lower limits, they allow us to perform comparisons in a self-consistent way. We also consider the r band alone as it is often used as a proxy for bolometric luminosity. The r band absolute magnitude is calculated as M_λ = m_λ − μ − A_λ − K_corr where μ is the distance modulus, A_λ is the Milky Way extinction in the corresponding wavelength and K_corr is the K-correction approximation K_corr = −2.5 × log(1 + z). We do not consider host extinction for any of the events as they have been deemed negligible or has not been considered by the different authors analyzing these events. With the aforementioned procedure, the r band light curve and the pseudo-bolometric luminosity light curves are presented in Fig. 1 and the peak value in Table 1. We can see that SN 2017egm is the most luminous event in the considered sample, closely followed by SN 2020wnt and SN 2018bsz.

Fig. 1. Comparison of the r band absolute magnitude (left) and the pseudo-bolometric (using the gri bands; right) luminosity light curves for the objects of our samples. The light curves have been aligned to r band peak; see Sect. 2.2 for more details. Fewer time bins appear in the luminosity panel because of the unfulfilled requirement of simultaneous g − r − i observations in some bins.

Given the diversity on the different light curves shape and spectral features, individual analysis conclude that most events should have at least some degree of CSM interaction. This does not exclude the possibility of a central engine in some cases. Gamma-ray emission in the giga-electronvolt to tera-electronvolt range is predicted for both the CSM interaction powered and magnetar powered scenarios, although no clear association has been published (see Renault-Tinacci et al. 2018, and references therein) until the recent case of SN 2017egm that is discussed in this work.

2.3. Fermi-LAT data analysis

The Fermi-LAT is a pair conversion γ-ray instrument operating in the 30 MeV to higher than 500 GeV energy range surveying the full sky since 2008 (Atwood et al. 2009), which makes it the ideal high-energy SN follow-up telescope. To study our sample of SLSNe we performed an analysis on a fixed window of one year after each SN discovery to minimize the number of trials. Then to characterize the emission we performed two different types of Fermi-LAT light curve analyses in the 100 MeV–100 GeV energy range. A first one spanning the entire dataset available (16.5 years; 2008-08-05 to 2025-02-05) with a time bin of 6 months using Pass8 R3 SOURCE events (Atwood et al. 2013; Bruel et al. 2018). Then a second light curve with finer 15-day time bins over 20 months around the SN (2 months before and 18 months after) is produced for the interesting candidates. As the LAT point-spread function (PSF) is highly energy-dependent, varying from several degrees width in the 100 MeV to 1 GeV band down to ∼0.1° for E > 10 GeV (Abdollahi et al. 2020), this tends to dilute faint signals and reduce the source significance in the low-energy band. To mitigate this issue, we used the summed likelihood method and jointly fit events with different angular reconstruction quality (the so called PSF0 to PSF3 event types⁷). We performed a summed likelihood with 4 components: three components (PSF1, PSF2, and PSF3) in the 100 MeV–1 GeV energy range, with a zenith angle cut of < 90° and a single component above 1 GeV with all event types and a broader maximum zenith angle cut of 105°. We performed a binned analysis within a region of interest (ROI) of 10° ×10° centered at the SLSN optical position with 8 energy bins per decade, and spatial bins of 0.05°. The data reduction was performed using the LAT Fermitools version 2.4.0⁸ and fermipy (Wood et al. 2017) version 1.4. To take into account the energy dispersion in the likelihood analysis⁹, we used the science tools parameter edisp_bins = −2. This means that two energy bins above and below the analysis energy range will be added when evaluating the model. The Galactic diffuse emission was modeled by the standard file gll_iem_v07.fits (Acero et al. 2016) and the isotropic diffuse emission were described by the tabulated model in iso_P8R3_SOURCE_V3_v1.txt. The models are available from the Fermi Science Support Center (FSSC)¹⁰. This setup is common for all all three analysis presented in this study.

To evaluate the significance of a putative γ-ray counterpart to the SLSN, we added a point source to the sky model at the optical SN position with a simple power-law spectrum E^−Γ with a fixed spectral index of Γ = 2. The statistical significance of the signal is estimated using TS = 2(lnℒ₁ − lnℒ₀), where ℒ₀ and ℒ₁ are the likelihoods of the background (including catalog sources, null hypothesis) and the hypothesis being tested (source plus background; see Mattox et al. 1996).

To account for background sources within the ROI our starting point is the latest release of the Fermi-LAT catalog (4FGL-DR4), which is based on 14 years of data (Abdollahi et al. 2022; Ballet et al. 2023). We added all sources up to a distance of 15° from the center of the ROI. We then used the optimize function of fermipy, which iteratively optimizes the parameters of the ROI in a multi-step process¹¹. Finally we performed another fit leaving only the normalization of the nearby sources (< 3°) free to vary, together with the isotropic and diffuse components. To check for potential new sources, we computed a TS map over the entire dataset period (16 years). We found no significant excess (TS > 25) in these residual TS maps.

Once our best ROI model was obtained, we performed our light curve analysis with the lightcurve¹² function of fermipy with 6-month and 15-day time bins where the Galactic diffuse background normalization is frozen to the best-fit value of the full time period. In order to accommodate variable sources in the ROI over time (other than our target), the normalization was set free for any source reaching TS > 16 in any time bin. Because most time bins yield TS ∼ 0, we computed upper limits using a Bayesian method following the approach of (Helene 1983). In particular, we employed the calc_int method of the IntegralUpperLimit class provided in the Fermi Science Tools. Hereafter, upper limits at a 95% confidence level are reported for detections with TS < 4 and flux points otherwise.

[13]Note that the bolometric luminosity used here is slightly different from the one shown in Fig. 1 because it includes more bands. We do this because here we want to characterize its energy, while in the other case we want a fair comparison with other events that have fewer observed bands.

3. Results

To minimize the number of trials, we first tested a simple Fermi-LAT analysis with the aforementioned setup in which we computed the significance of a γ-ray source at the SN optical position, in a fixed time window (T_SN ⇒ T_SN + 12 months), and with a fixed spectral index of 2.0. The time window is motivated by the fact that the magnetar model predicts the γ-ray emission to peak in the first 12 months post-explosion (Vurm & Metzger 2021) and that dense CSM shells close to the star are needed to explain the early part of the superluminous light curve in the optical (e.g., Lin et al. 2023). The fixed spectral index is motivated by the fact that the standard diffusive shock particle acceleration (DSA) mechanism produces a particle population with an index close to 2 and that in the magnetar model the expected SED is rather flat (i.e., also with Γ ∼ 2) in the first months (Vurm & Metzger 2021). The results of this first analysis are reported in Table 1. We note that only SN 2017egm shows a significant γ-ray emission within the first year with a TS of 25.7 corresponding to ∼5σ pretrial for one degree of freedom (the normalization) and 4.7σ post-trial considering that six SLSNe were tested. When the spectral index is let free (2 degrees of freedom) for SN 2017egm, the best fitted index is 2.1 ± 0.2 for a TS of 28.1, which also corresponds to ∼5σ.

For SN 2019ieh, only a small hint of a signal is observed and is discussed further in the following section with a 6-month time bin light curve. To evaluate the impact of the fixed spectral index hypothesis, we changed it to 1.5 (2.5), which resulted in a TS of 9.6 (4.3) for SN 2019ieh.

3.1. Long-term light curve of the SLSN sample

While the γ-ray emission in the magnetar model is expected in the first year after the SN, we also considered a time range covering a larger dataset from the Fermi launch in August 2008 to August 2024 and a time binning of 6 months. This provides the historical light curve of the SN host galaxy. For example, faint signals over the 16 year period could indicate possible active galactic nucleus (AGN) activity in the host galaxy. On the other hand, the detection of gamma-ray flux within the year after the SN was observed optically, combined with the lack of historical variability of the host galaxy, would strengthen the possibility that, at least temporally, the γ-ray signal could be related to the SN.

The light curves resulting from our analysis are presented in Fig. 2. One can see that besides SN 2017egm, no time bin with TS > 4 is observed for other targets, which include SN 2018bsz at a distance of 111 Mpc, the nearest SLSN I recorded to date (Anderson et al. 2018), SN 2020wnt, which shows the second-most luminous light curve (see Tab. 1), and SN 2021adxl, the closest event in the SLSN II sample of Pessi et al. (2025). For SN 2019ieh, the light curve with 6 month time bins confirms that there is no significant emission. Given the limited statistics available, the current data cannot definitively determine whether the observed signal stems from a genuine astrophysical source or is merely attributable to random statistical fluctuations. In addition, given the fact that SN 2019ieh is the least luminous of our sample we do not consider it further in this study. A discussion of the possible reasons for the lack of γ-ray detection for these targets is presented in Sect. 6.

Fig. 2. Luminosity light curves in the 100 MeV – 100 GeV energy range over 16 yrs for each SN of our sample from the Fermi launch to August 2024 with a time bin of 6 months. A flux point is shown when the TS > 4, otherwise upper limits at the 95% confidence level are reported. In all time bins, the spectral index of the tested source is fixed to 2. The SN discovery date is also indicated. For a comparison across the sample, the derived Fermi-LAT flux is transformed to luminosity using the distance indicated in Table 1.

3.2. Refined light curve of SN 2017egm

Given the light curves presented in the previous section, SN 2017egm is the only target for which a γ-ray signal is observed and it will be the main focus of the rest of this work. To get a more detailed understanding of the time evolution of the γ-ray signal, we performed a refined light curve Fermi-LAT analysis focused around the SN discovery time. The time binning is refined to 15-day time bins starting 2 months before the SN explosion and extending to 18 months after. For this finer-binned analysis, we used the same analysis configuration (summed likelihood with different PSF types) as presented in Sect. 2.3. This more finely binned light curve confirms the detection of a delayed γ-ray emission that peaks approximately 120 days after the SN explosion, as is shown in Fig. 3.

Fig. 3. Left: comparison of the optical bolometric luminosity13 taken from Lin et al. (2023) and the γ-ray 100 MeV–100 GeV luminosity (this work) for SN 2017egm along with the TS in each 15-day time bin and with larger time blocks presented in Sect. 3.2 in black. Middle and right: γ-ray and optical luminosity observations and predictions from the magnetar model tailored to SN 2017egm by Vurm & Metzger (2021) for different nebula magnetization parameters (see Sect. 5.1). Hereafter the T0 value is taken at the discovery date 23 of March 2017.

In the light curve shown in Fig. 3 (left panel) no signal is detected within the first ∼50 days after the explosion. Then the γ-ray emission rises, reaching a peak significance of TS ∼ 22 at T ∼ 130 days post-explosion. The flux is consistent with the upper limits derived by Acharyya et al. (2023), whose analysis differs in the exact time range employed and the optimization of the data selection.

To assess the duration of the peak and possible substructure, we applied the Bayesian block algorithm (Scargle et al. 2013), following the likelihood formalism developed by Kerr (2019). Within this framework, the ROI model was employed to compute the probability for each photon in the ROI to be associated with the new source with gtsrcprob¹³, and used as a weight to compute the number of blocks in our light curve (N_b). The addition of further blocks was penalized with a prior $\propto N_{\mathrm{b}}^{\gamma }$, where γ parametrizes the false-positive rate. We computed this rate by redistributing photons randomly among blocks, and found that γ = 5 approximately corresponds to 2σ. A total of three blocks are found: two non-detections and a γ-ray signal well represented by a single block from 57939 MJD (2017-07-05) to 58051 MJD (2017-10-25) for a total duration of 112 days. No significant substructure is found for a variability threshold set at 2σ. We computed upper limits for the blocks with non-detections. The first corresponds to a single bin prior to the γ-ray emission (from MJD 57835 to 57939), while the last block is divided in two bins (MJD 58051 to 58192 and MJD 58239 to 58482) as per the coverage limitations owing to recovery from the failure of one of the solar array rotator drives on Fermi¹⁴. The resulting upper limits are shown in Fig. 3.

For illustrative purposes, the duration of the Bayesian block is also shown in Fig. 4 on top of the photon distribution as a function of time and angular distance to the SN optical position. The probability of being associated with the SN and the reconstructed energy for each photon is also added. As was expected, the density of near-target photons with high association probability is higher within the Bayesian block relative to the full analysis period.

Fig. 4. Photon angular distance to SN 2017egm optical position as a function of time. Using the best-fit model of the ROI, the probability that a given photon is associated with the SN is illustrated by colors, while the size indicates the photon energy. The γ-ray time interval defined by the Bayesian block algorithm is illustrated with the shaded area. Coverage limitations due to the failure of one of the solar array rotator drives on Fermi resulted in the loss of approximately 18 days of data in March and April 2018.

3.3. Spectral analysis and source localization

Using the time block defined earlier, we carried out a spectral study using 112 days of data with the same Fermi-LAT analysis configuration (see Sect. 2.3). Over this short period of time, and at an extragalactic position, the level of Galactic diffuse emission cannot be well constrained, and its normalization is fixed to the best-fit value derived from the 20-month analysis (norm = 0.89). The isotropic diffuse emission and the normalization of the brightest source 4FGL J1015.0+4926 in the ROI was set free. The parameters of other sources were fixed to the fit obtained over the 20-month period. Assuming a power-law spectral model, the likelihood analysis results in a significance of TS = 33 for SN 2017egm when the position is fixed at the SN optical position. With two degrees of freedom for the source (normalization and spectral index), this corresponds to a 5.4σ detection.

We note that in the Bayesian block time interval, one photon of high energy (E = 6.8 GeV) has a high association probability (0.94) as shown in Fig. 4 and Table A.1. To test the robustness of our detection against a single photon, we performed a likelihood analysis where this single photon is removed. This resulted in a TS = 25 (instead of 33) indicating that most of the significance does not depend on a single photon.

When fitting the source localization using the fermipylocalize command, the likelihood improvement is marginal (Δlogℒ  = 0.5). This indicates that the γ-ray signal is statistically compatible with the optical position. The best localization is at RA_J2000 = 154.81°, Dec_J2000 = 46.41° (∼0.05° from the optical position; see also Sect. 4) with a 95% error radius of 0.12°. This localization is compatible with the position reported in Li et al. (2026) with a slightly smaller 95% error compared to their study (0.18°).

The residual TS map (where only SN 2017egm is excluded from the model) is presented in Fig. 5 (top panel) and shows no significant residuals in the 10° ×10° region around the source. The SN best-fit spectral model parameters are Γ = 2.17 ± 0.23 and a normalization N₀ = (0.82 ± 0.22)×10⁻¹² cm⁻² s⁻¹ MeV⁻¹ for a reference energy of 1 GeV. We tested for spectral curvature with an exponential cutoff power-law model, which resulted in a marginal likelihood improvement of Δlogℒ ∼ 3. The SED shown in Fig. 5 (bottom panel) was estimated by computing the photon flux in each energy interval, assuming a power-law shape with a fixed photon index of Γ = 2 for the source. A 95% confidence level upper limit was computed when TS < 1 in a given energy bin. Spectral data points are given in Table A.2.

Fig. 5. Top: residual TS map in the 0.1–100 GeV energy range in the time interval defined by the Bayesian block algorithm (2017-07-05 to 2017-10-25). The plus symbols represent the sources from the 4FGL-DR4 catalog and the SN optical position. Bottom: γ-ray luminosity SED of SN 2017egm in the same time interval. The spectral predictions from the magnetar model are shown for different nebula magnetization parameter (see Sect. 5.1 for more details).

3.4. Comparison with past studies

Our detailed study confirms the evidence found by Li et al. (2026). We note that there are some variations in the analysis setup such as energy range, time range, and spatial binning size that could explain some of the differences in TS that are observed. In particular their 1-year analysis (shown in their Fig. S2) starts on August 4, 2017, and yields a TS = 23.8. As our Bayesian block analysis yields a starting date on July 5, 2017, we can assume that some signal is lost by starting in August. This difference and the fact that we carry out a joint likelihood with PSF event type (i.e., a more sensitive analysis) might explain the slight differences observed. Concerning the spectral parameters obtained in this study (Γ = 2.17 ± 0.23), they indicate a marginally harder spectrum compared to those reported by Li et al. (2026) of Γ = 2.6 ± 0.4. This difference may be attributed to the differing energy ranges employed; 0.5–500 GeV in the aforementioned study, whereas our analysis was conducted within the 0.1–100 GeV range. If we reproduce the same energy range of their analysis, we obtain similar results (Γ = 2.6 ± 0.3). Regarding the analysis by Acharyya et al. (2023), similar variations could explain the observed differences. In particular for faint sources, the coarse spatial bin size (0.1°) can significantly degrade the source significance as illustrated in Table 2 of Acero et al. (2022) in the case of the faint γ-ray signal of the Kepler supernova remnant. We also note that Li et al. (2026) investigated the specific setup of Acharyya et al. (2023) and were able to reproduce a similar TS.

4. Exploring other multiwavelength associations

While the giga-electronvolt association with SN 2017egm is certainly enticing, the 0.12° uncertainty in the localization (95% confidence level) requires further evaluation. Most extragalactic giga-electronvolt transients are identified as flaring blazars; therefore, AGN activity could be a suitable interpretation for the excess. Li et al. (2026) searched for optical flares in the Gaia archive (Gaia Collaboration 2023) and cross-matched nearby radio sources with an analysis of Chandra data, finding three alternative counterparts (see Table II in Li et al. 2026). One corresponds to a source within the host galaxy NGC 3191 (Radio Source 1, coincident with WR 187; see Brinchmann et al. 2008) and another to the galaxy MCG+08-19-017 (Radio Source 2). None of them show nuclear activity, and the third source (SDSS J101815.31+462626.6, Radio Source 3) falls outside our 95% containment region.

However, the region contains multiple, relatively faint X-ray sources reported in the XMM-Newton DR14 and Chandra point source catalogs CSC 2.1 (Webb et al. 2020; Evans et al. 2024) with no known radio counterpart (see Fig. 6). For example, Li et al. (2026) discuss the optical emission of two quasars (SDSS J101918.93+462301.4 and SDSS J101858.02+461812.3) that are spatially consistent with detected X-ray sources (4XMM J101918.8+462301 and 4XMM J101858.0+461812, respectively). These sources have no reported variability in the XMM-DR14 catalog. The brightest among all X-ray sources within our localization error is the quasar 2CXO J101839.5+462315 (SDSS J101839.48+462315.8, z = 0.941; Albareti et al. 2017), which is only reported to have mild variability in the CSC 2.1 catalog. Notably, the largest flux is found at ∼9.2 × 10⁻¹⁴ erg/cm²/s (ObsIDs 19034 and 19035; 2017-09-17 and 2017-11-09, respectively), compared with ∼5.6 × 10⁻¹⁴ erg/cm²/s derived in the only other existing observation (ObsID 19036; 2018-05-21). No optical variability is found in quasi-contemporaneous ATLAS data (Tonry et al. 2018), although we note that the limiting magnitude of its forced photometry is close to the brightness of the quasar. Only mild variability (Δm ≤ 1) is found by Gaia, PTF, or ZTF (Rau et al. 2009; Masci et al. 2019; Gaia Collaboration 2023) – considering observations prior or after the giga-electronvolt excess peak. Furthermore, it is worth noting that SDSS J101839.48+462315.8 lacks a known radio counterpart in any of the major radio surveys (VLASS, LoTSS, RACS, FIRST, NVSS, and SUMSS; see Flesch 2024). According to the Fermi-LAT AGN catalog by Ajello et al. (2020), the vast majority of γ-ray AGNs exhibit radio fluxes exceeding approximately 5 mJy (refer to their Fig. 14). Considering that current surveys achieve sensitivity levels in the millijansky range (and down to sub-millijansky levels for VLASS, Lacy et al. 2020), the apparent radio-quiet nature of SDSS J101839.48+462315.8 suggests that it does not conform to the typical characteristics of a Fermi-LAT blazar.

Fig. 6. Multiwavelength maps of the region around the giga-electronvolt excess: radio (VLA; left) and X-rays (XMM; right). The white marker represents the best-fit position of the giga-electronvolt excess, and the dotted curve is its localization error ellipse at a 95% confidence level. The cyan and green markers locate the nearby galaxies NGC 3191 and MCG+08-19-017, while the red and blue correspond to NVSS J101948+462332 and SDSS J101839.48+462315.8, respectively. White circles mark the quasars SDSS J101858.02+461812.3 and SDSS J101918.93+462301.4 discussed by Li et al. (2026).

The last relevant candidate might be the unidentified VLA source NVSS J101948+462332 (4.9 mJy at 1.4 GHz; Condon et al. 1998). But the source has no X-ray or optical counterpart, and no variability was found in radio (Ofek & Frail 2011). Although the radio source lies outside Chandra’s field of view, XMM observations set a limit at F(0.2 − 2 kev) < 3.9 × 10⁻¹⁵ erg/cm²/s (Obs_id: 08437401; König et al. 2022) – for N_H = 10²⁰ cm⁻² (HI4PI Collaboration 2016) and assuming Γ = 2. As in the previous case, NVSS J101948+462332 is an unlikely association given the low X-ray flux compared with the known population of giga-electronvolt-bright AGN (Ajello et al. 2020).

Overall, although the giga-electronvolt excess is spatially consistent with the position of several quasars, the existing multiwavelength monitoring does not strongly support any association other than SN 2017egm. At the same time, the poor coverage during the Bayesian block time interval prevents an unambiguous rejection of the associations discussed.

5. SN 2017egm γ-ray emission: CSM interaction and magnetar nebula scenarios

As the first SLSN to be detected at γ-ray energies, SN 2017egm represents a new window to shed light on the mechanism powering these objects. To understand the nature of the γ-ray emission, the temporal and spectral results from this SN are discussed in the magnetar and CSM interaction frameworks.

5.1. Magnetar-powered γ-ray emission

In the magnetar model, the rotational energy from the millisecond magnetar is converted into a wind of electrons and positrons that inflates a magnetar wind nebula of highly energetic particles. High-energy photons from the nebula will be injected into the ejecta via IC and synchrotron processes, where they will be partly absorbed and reprocessed into thermal optical/UV radiation. The absorption implies mechanisms such as the Bethe-Heitler process (photon-matter interaction) and γ − γ pair production on thermal (optical/UV) or nonthermal (X-ray) photons that lead to an attenuation of the γ-ray flux. These interactions govern the energy transfer and the escape of high-energy radiation. This mechanism, where high-energy photons from the nebula convey energy to the ejecta, is a fundamental aspect of the magnetar model for maintaining the high luminosity of SLSN over an extended period.

To accurately describe the thermalization process and the optical depth as a function of both energy and time in a self-consistent way is a challenging task due to the interplay of all the interaction processes and the evolving properties of the ejecta. To tackle this issue, Vurm & Metzger (2021) employed a three-dimensional Monte Carlo time-dependent numerical radiative transfer model that tracks the production, transport, and interactions of photons and electron-positron pairs throughout the nebula and the ejecta. Of particular interest for this study, the predicted observables from their model are the optical luminosity, the inferred escaping γ-ray luminosity and associated spectral distribution.

To adapt their model to the specific case of SN 2017egm, the system parameters (magnetic field strength of 1 × 10¹⁴ G, pulsar period at birth of 5 ms, and ejecta mass of 3 M_⊙) have been fixed from the optical light curve from Nicholl et al. (2017). Once these parameters are set, the remaining main free parameter in their model is the nebula magnetization parameter, ε_B. This is the ratio of the magnetic energy over the injected magnetar energy (see Eq. (14) in. Vurm & Metzger 2021). In a more magnetized nebula, a greater fraction of the energy is radiated via synchrotron emission. This emission can be more readily absorbed and thermalized by the ejecta than the higher-energy IC emission, which is more prone to escape without thermalizing (see Sect. 2.2 and Eq. (17) of Vurm & Metzger 2021). As can be seen in the models shown in Fig. 3 (middle and right panels), the direct consequence of this is that for a nebula with high magnetization, a greater fraction of the engine power is thermalized producing greater late-time optical luminosity and lower levels of γ-rays.

From the optical and giga-electronvolt light curves, and the γ-ray spectrum, a magnetar model with low nebular magnetization is plausible with ε_B ∼ 10⁻⁶ or smaller. This is in particular required not to largely overshoot the optical light curve at t > 250 days. As is discussed in Vurm & Metzger (2021), such a low magnetization value (compared to > 10⁻² for the Crab nebula but at a very different age) implies the existence of an extremely efficient magnetic dissipation mechanism in the nebula. An interesting alternative scenario to reduce the late-time optical luminosity is that the spin-down luminosity decays faster than the typical magnetic dipole rate ∝t⁻²; for example, due to a growing magnetic dipole moment or fall-back material accretion (Metzger et al. 2018).

In the magnetar scenario, given a sufficiently bright source, the γ-ray observations at giga-electronvolt at late-time would allow one to trace the time evolution of the spin-down luminosity by measuring the flux decay after the light curve peak. Unfortunately this is not possible for SN 2017egm, as the upper limits derived at later times are not sufficiently constraining. The role of tera-electronvolt observations will be discussed later in Sect. 7.

5.2. A hadronic CSM interaction scenario

As an alternative to the magnetar scenario, it has been proposed that the optical component and its late time bumps can be interpreted as ejecta interactions with a multi-component CSM (Zhu et al. 2023; Lin et al. 2023). In this scenario, hadrons from the engulfed CSM would be accelerated via diffusive shock acceleration at the forward shock, interact with the surrounding CSM and produce γ-rays through hadronic channels – predominantly π⁰-decay. In the following, we do not consider contributions from the reverse shock and will disregard leptonic contributions – presumed to be subdominant given the high CSM densities.

Based on the CSM characterization of Lin et al. (2023) in four distinct shells (Appendix B), we constructed a density profile, ρ(t), encountered by the shock assuming an expansion parameter m = (n − 3)/(n − s) = 0.75 (n = 12, s = 0) and an initial velocity of V_s, 0 = 9.8 × 10³ km/s (see Eq. (B2) in Chatzopoulos & Wheeler 2012; Lin et al. 2023). Instead of a stellar wind profile, a constant density is assumed within each shell – a scenario also considered by Wheeler et al. (2017) for the optical peak. Diffusive shock acceleration would occur at the forward shock, although it may be suppressed if the shock is radiation dominated (see, e.g., Ito et al. 2018; Levinson & Nakar 2020). Here two scenarios are considered: (1) that particle acceleration starts at least one day after the explosion (t_i = T₀ + 1 d; as in Tatischeff 2009), or (2) hadrons can only gain energy beyond R_s(t_i)∼R_br, after shock breakout (a conservative assumption; see, e.g., Murase et al. 2011; Katz et al. 2012; Giacinti & Bell 2015). The energy-conversion efficiency of accelerating cosmic rays can be parametrized as η in the cumulative energy, E_CR (i.e., as a fraction of the shock’s kinetic energy), which is estimated as

$$\begin{array}{c}\hfill E_{\mathrm{CR}}=\eta {\displaystyle {\int }_{t_i}^{t_f}}\frac{1}{2}\rho (t)V_{\mathrm{s}}{(t)}^{3}4\pi R_{\mathrm{s}}{(t)}^{2}dt,\end{array}$$(1)

where R_s(t) is the shock radius and $V_{\mathrm{s}}=d{R}_{\mathrm{s}}/dt=V_{\mathrm{s},0}{\left[\frac{t}{1\mathrm{day}}\right]}^{m-1}$ (see, e.g., Tatischeff 2009; Martí-Devesa et al. 2024). Shock breakout from the first shell occurs when the optical depth is τ_opt ∼ c/V_s(t_br), resulting in a breakout time of t_br ∼ 14 d, assuming an opacity of κ = 0.2 cm²/g (Lin et al. 2023).

The energy budget, E_CR, should be compared with the total energy stored in protons required to explain the giga-electronvolt excess (typically expressed as W_p). We assume that during the time interval defined previously through Bayesian blocks (see Sect. 3.2) the average proton distribution is well reproduced by a power law with an exponential cutoff in momentum space, i.e.,

$$\begin{array}{c}\hfill \frac{d{N}_{\mathrm{p}}}{\mathit{dE}}=\beta N_0{\left(\frac{E}{E_0}\right)}^{-p}exp-\left(\frac{E}{E_{\mathrm{cutoff}}}\right),\end{array}$$(2)

where N₀ is a normalization factor, E₀ the scale energy, E_cutoff the cutoff energy, p the proton index, and β the particle velocity, v/c – correcting the spectral shape of a proton population below ∼1 GeV when described in kinetic energy (Dermer 2012). At detection time, the maximum energy is not limited by proton-proton collisions. Given the lack of curvature or simultaneous VHE observations (Acharyya et al. 2023), we consider again two extreme scenarios for E_cutoff – at 1 TeV and 1 PeV. We employ the radiative process and fitting package Naima to perform an MCMC fit of a π⁰-decay component (Kafexhiu et al. 2014; Zabalza 2015) to the Fermi-LAT spectrum. From the previously derived ρ(t) profile, we calculated the average density encountered by the shock during the Bayesian block time interval, which is 1.5 × 10⁻¹⁴ g cm⁻³. The fit results in a proton distribution with $p=2.5_{-0.3}^{+0.4}$, with a total energy content of $W_{\mathrm{p}}(>0.1\mathrm{GeV})=3.8_{-1.7}^{+8}\times {10}^{49}$ erg ($\eta \sim 0.{05}_{-0.02}^{+0.09}$). While energetic considerations are satisfied (Fig. 7, left panel), it is apparent that the γ-ray emission should trace the target medium in a CSM scenario. With the cosmic-ray efficiency and spectral index derived above (η = 0.05, p = 2.5), we can estimate the γ-ray flux that would be produced by the succession of CSM shell interactions and produce the corresponding γ-ray light curve. Absorption processes (Bethe-Heitler or pair-production) are not considered here, but these do not affect our conclusions and may only impact the flux at earlier times for the LAT energy range (Vurm & Metzger 2021). The resulting time-evolving γ-ray emission is shown in Fig. 7 (right panel) and we note that the narrow CSM shell structures provide a γ-ray emission time stamp not matched by observations (a time delay between the predicted γ-ray emission from a given shell and the observed γ-ray signal). As this CSM model was built by Lin et al. (2023) to reproduce the features of the optical light curve, our observed mismatch can be interpreted in the CSM model as a temporal offset between the optical and the γ-rays.

Fig. 7. Left: energy conversion into cosmic rays as a function of time for different η values (Eq. (1)), and compared with the best fit obtained from Naima. Vertical lines represent the crossing times of the inner radius of each shell. Scenarios for ti = 1 d (green) and ti = tbr (black) are displayed, both assuming Ecutoff = 1 PeV. Right: γ-ray flux expected from the interaction of each shell (η = 0.05) compared with the Fermi-LAT flux (100 MeV–100 GeV). The arrow represents the Bayesian block interval of the emission peak, while the TS is shown as shaded green bars. Both panels display the relevant time interval for the CSM shell model.

We note that in dense CSM shocks, the optical and γ-ray emissions originate from fundamentally different physical processes. Optical radiation is produced by the shock-powered thermal X-rays that are efficiently absorbed and reprocessed, while γ-rays arise from the nonthermal proton–proton interactions at the forward shock. In principle, this difference in emission sites creates the possibility that the optical and γ-ray light curves may not track each other in time, because the photospheric radius and the forward-shock radius can decouple if the diffusion time becomes long enough. However, unless the diffusion time greatly exceeds the dynamical time (an optical depth of τ ≫ c/V_s), the reprocessed luminosity still emerges broadly at the same moment with the nonthermal radiation from the accelerated particle population. Such conditions would imply a radiation-mediated shock (τ ≫ 1), which cannot efficiently accelerate cosmic rays, as was explained earlier. Therefore, the regime that would allow long optical delays is incompatible with the one accelerating particle in a collisionless shock and producing giga-electronvolt photons.

In addition to the timing mismatch argument, one can compare the energy radiated in optical and in γ-rays at ∼100 days post-SN explosion, which in our case yields L_γ/L_opt ∼ 1 for simultaneous observations (see Fig. 3). As the shock converts > 90% of its kinetic energy into thermal energy and < 10% into cosmic rays, and then only a fraction become γ-rays, one would expect at most L_γ/L_opt ∼ 10⁻² − 10⁻¹, which is one or two orders of magnitude lower than the measured ratio. In other words, as the thermal channel is far more efficient than the nonthermal one, one should expect a higher optical luminosity than the γ-ray one. This ratio can be contextualized with systems in which both the optical and γ-ray outputs are known to be powered by CSM shock interactions. For the nearby Type II SN, SN 2023ixf–one of the closest core-collapse events in the past decade–Fermi-LAT non-detections constrain the γ-ray–to–optical luminosity ratio to (L_γ/L_opt < 10⁻²) at epochs somewhat earlier than those considered here (see Fig. 9 of Martí-Devesa et al. 2024). Likewise, in the recurrent nova RS Oph, where CSM interaction is firmly established, the Fermi-LAT detections reveal a nearly constant ratio (L_γ/L_opt ≃ 2 − 3 × 10⁻³) over the first ∼50 days (Cheung et al. 2022). These values are representative of shocks in dense environments and highlight the strong dominance of the thermal (optical) channel over the nonthermal (γ-ray) in interaction-powered transients. As a conclusion, the CSM model faces two fundamental challenges: 1) the time delay cannot be simply explained by diffusion effects, and 2) the CSM model should produce much more optical emission than γ-rays due to the higher efficiency of the thermal channel.

6. Discussion

6.1. The specific properties of SN 2017egm

An open question remains regarding why only SN 2017egm is detected in our sample and not the other sources at similar distances or closer. Among the SLSN I, SN 2017egm is the only source in our sample exhibiting “W-shape” O II features and He features (Zhu et al. 2023), while SN 2018bsz and SN 2020wnt stand out for their C II features (we do not discuss SN 2109ieh as it is not as luminous and this could be the sole reason why γ-rays are not detected). The presence of this O II feature is thought to be a consequence of excitation of CNO layers by nonthermal radiation for example from the magnetar wind nebula (Mazzali et al. 2016). These spectroscopic differences point toward different progenitor characteristics (ejecta mass, presence of a central engine, differences in the rotational velocity of the magnetar, and so on) that might impact the time evolution of the opacity and the expected γ-ray flux. Whether SN 2017egm is representative of SLSNe in general is not yet clear and a larger SLSN population would be needed to further investigate this question.

SN 2020wnt, the second-most luminous event in our sample, has been associated both with a radioactive nickel scenario (Gutiérrez et al. 2022) and a magnetar scenario (Tinyanont et al. 2023), with the latter providing a better fit to the observed light curve, and thus being preferred. Given the spectral evolution, Tinyanont et al. (2023) suggests that the magnetar is hidden inside the optically thick ejecta near peak, and only exhibited magnetar spectroscopic signatures at later epochs (∼1 year after peak). This could be due to a possibly larger ejecta mass that might increase the γ-ray opacity and delay the time window when γ-rays could leak from the ejecta to a moment when the spin-down luminosity had already significantly decreased and is no longer detectable.

6.2. Limitations of the CSM and magnetar models

All the events in our sample show some degree of interaction (light curve bumps and/or narrow spectral lines). While some levels of γ-ray emission can be expected in the CSM-shock interaction, we note that it is not guaranteed that the conditions for efficient particle acceleration (requiring a collisionless shock) are reached in the early times. First, the shock will be radiation-mediated, and only then, when photons escape the system, a collisionless shock can form and particle acceleration becomes possible (Murase et al. 2014). In addition, in a very dense CSM and photon radiation field, the γ − γ absorption process can play an important role and drastically decrease the expected γ-ray radiation at E > 10 GeV (see, e.g., Cristofari et al. 2022).

The models presented in Fig. 3 show good agreement with the initial light curve decline and the γ-ray light curve but fail to capture the bumpy late-time tail of the optical emission. Thus, a possible solution consists of a hybrid model in which the early time light curve is dominated by the magnetar model and the late time, in particular the bumpy features, by CSM interaction. This scenario requires a very low magnetization of the nebula (ε_B < 10⁻⁶) or a spin-down luminosity decaying faster than a magnetic dipole rate of t⁻², which would need to be explored in more depth theoretically, in a scenario where infalling accreting material could play a role (Metzger et al. 2018).

Alternatively, it has been recently proposed that the bumps in the optical light curve could be related to modulation of the central engine luminosity in particular due to an infalling accretion disk undergoing Lense-Thirring precession as introduced in Lense (1918)– see also the review by (Mashhoon et al. 1984). The general concept is that around the newly formed magnetar a small, tilted accretion disk from failed ejected matter is formed. A misalignment between the accretion disk axis and the magnetar spin axis induces a Lense-Thirring torque that causes the disk to precess. This modulates the magnetar emission being injected into the ejecta, and therefore the amount of optical emission (see Farah et al. 2025). Note that in this model the magnetar model explains the bulk of the luminosity and the precession of the accretion disk the small modulation. This mechanism was successfully applied to model the light curve bump structures observed in optical in the SLSN 2024afav (Farah et al. 2025). Whether a similar scenario can quantitatively be applied to SN 2017egm remains to be seen, but such models would have the benefit of providing a self-consistent picture.

7. Perspectives for tera-electronvolt detection

As the γ − γ optical depth in the pulsar wind and ejecta decreases with time, the γ-ray window opens to > 10 GeV energies within the first year after the SN and the peak energy of the escaping γ-rays shift to higher energies as can be seen in Fig. 8 (top panel). However, the optical depth still remains high above tera-electronvolt energies (see Fig. 4 of Vurm & Metzger 2021) meaning that the window for observations by atmospheric Cherenkov telescopes is mostly restricted to the ∼10 GeV–1 TeV energy band.

Fig. 8. Top: Comparison of CSM and magnetar (magnetization εB = 10−6) spectral models of SN 2017egm at a distance of 135 Mpc together with Fermi-LAT (this work) and VERITAS (Acharyya et al. 2023) observations. In the CSM model, the predicted flux is shown with or without γ − γ absorption (see Appendix C), highlighting the fact that no emission is expected above 50 GeV. A simulation of 50 hours of CTAO-North observations between 50 GeV and 10 TeV is also shown for the magnetar model (see Sect. 7 for more details). For clarity only detection flux points (where TS > 4) are shown and upper limits are omitted. Using the CTAO simulation shown above, below we display the significance of a SN 2017egm-like source at different distances in the magnetar model only.

An upper limit on the tera-electronvolt emission from SN 2017egm was reported by the VERITAS collaboration and is presented alongside the Fermi-LAT flux measurements in Fig. 8 (top panel). The VERITAS observations correspond to a livetime of 8.7 hours and a lower-energy threshold of 350 GeV. Regardless of whether the Fermi-LAT emission originates from a CSM shock interaction or is powered by a magnetar central engine, the predicted tera-electronvolt fluxes from both models lie well below the reported VERITAS upper limit. Notably, the magnetar model exhibits a spectral peak in the 100–200 GeV range, followed by a sharp decline at higher energies, making the VERITAS low-energy threshold of 350 GeV insufficient for detection in this scenario. However, given the better sensitivity in the sub-tera-electronvolt energy range of the Cherenkov Telescope Array Observatory (CTAO), we explore its sensitivity to detect SN 2017egm-like objects in the future in the magnetar model at different distances.

The CSM model is not simulated as its tera-electronvolt spectrum is severely hampered by the γ − γ absorption (see Appendix C) and no emission is expected above 50 GeV. For the magnetar case, as its model predicts a curved spectrum peaking at 100 GeV, we only simulated observations with CTAO North, which with four large-size telescopes (LSTs) is better optimized for the low-energy range. For the simulations of the expected significance, we used the package Gammapy (Donath et al. 2023) with version 1.3 (Acero et al. 2025). To simulate the magnetar wind-powered spectrum, we used the model of Vurm & Metzger (2021) for the magnetar emission at T = 300 days (ε_B = 10⁻⁶) as at T = 100 days no emission above 50 GeV is expected.

The spectrum is then loaded in Gammapy using the TemplateSpectralModel class. A 3D (RA, Dec, Energy) Gammapy dataset was constructed by assigning a point source with the aforementioned spectral model at the center of a 3° ×3° region of interest with a bin size of 0.02° and 15 logarithmically spaced energy bins from 50 GeV to 10 TeV. The instrumental background and responses were derived from the CTAO-North instrument response functions (IRFs¹⁵) : Prod5-North-20deg-AverageAz-4LSTs09MSTs.180000s.

Assuming 50 hours of observations, the resulting simulated flux points are shown in Fig. 8 (top panel) for the magnetar models, showing that a detection and a reasonable characterization of the spectral properties can be obtained with CTAO for such objects. Interestingly, in the CSM model as the γ − γ absorption is likely to kill any signal in the 50 GeV to 50 TeV energy range, a detection at tera-electronvolt energy of a SLSNe would be a strong argument in favor of the magnetar model. However, it is important to keep in mind that these predictions are model dependent and based on the sole case of SN 2017egm, which may not be representative of the population of Type I SLSNe.

The detection significance was estimated using an Asimov dataset (a data cube where each bin is set to the predicted value N_pred; see Cowan et al. 2011, for a formal mathematical description), which provides the expected median significance of detection without running intensive Monte-Carlo simulations. The significance is then the likelihood ratio of the 3D Gammapy analysis with or without the putative point source assuming a power-law model with a fixed spectral index of 2.5 and leaving only the normalization of the source free.

For the SN 2017egm magnetar model at a distance of 135 Mpc, a significance of TS = 30 (∼5.5σ) indicates that SN 2017egm could have been detected by CTAO with 50h of observations. To test the horizon of detectability we considered candidates in a distance range from 70 Mpc to 200 Mpc as shown in Fig. 8 (bottom panel). Our conclusion is that any SN 2017egm-like object in the magnetar model within 100 Mpc would yield a strong detection (> 10σ), while distances greater than 100 Mpc would exhibit fainter detections and be more subject to the choice of model parameters. Anyhow, this simulation shows that nearby SLSNe could represent a new class of source for CTAO. Future detections in the CTAO energy range could provide constraints on the central engine powering SLSNe, the magnetar properties, and, in particular, the magnetization of the nebula.

8. Conclusions

Based on a curated list of nearby Type I (hydrogen poor) and Type II (hydrogen rich) SLSNe, we performed a systematic search for γ-ray emission within the last 16 years of Fermi-LAT observations using a summed likelihood analysis employing the LAT PSF types, optimized to search for faint emission at the optical positions of the SLSN. Then we characterized the light curve shape with a refined analysis with time bins of 15 days and a Bayesian block algorithm for the case of SN 2017egm. Our main conclusions can be summarized as:

1. From the targets in our sample, only SN 2017egm is detected, confirming the report by Li et al. (2026). We found the likelihood TS to be TS = 26 in a fixed one-year time window post-explosion or TS = 33 in the Bayesian time block. Both TS are reported assuming a position fixed to the SN optical one. The optimal time window starts ∼50 days after the explosion and lasts for 112 days up to late October of the same year. Within that time frame, the γ-ray emission can be modeled with a power-law of index Γ = 2.17 ± 0.23 with no evidence of spectral curvature. The localization of the γ-ray source is statistically compatible with the SN optical position.

2. Our spectral and light curve results are visually compared against CSM and magnetar models. In the latter model of Vurm & Metzger (2021), the γ-ray emission from the magnetar wind nebula matches both the observed flux level and the light curve shape. In particular, the model predicts γ-rays to be visible only after ∼60 days, once the opacity of the ejecta has decreased, which is compatible with our findings. However, in order to explain both the optical and the γ-ray light curves, an extremely low magnetization of the nebula is required (ε_B < 10⁻⁶) or alternatively a spin-down luminosity decaying faster than a magnetic dipole rate.

3. Concerning the CSM model, assuming the multiple CSM shell model of Lin et al. (2023), we estimated the expected γ-ray flux stemming from each shell in a hadronic scenario. While the flux level is compatible with energetic considerations, the timing of the observed γ-ray signal and the ratio L_γ/L_opt ∼ 1 are inconsistent with theoretical expectations. This ratio is also not in line with measurements in other CSM-dominated objects such as novae and SNe, which yield L_γ/L_opt < 10⁻².

4. Based on the previous models, we conclude that the γ-ray emission is well described by the magnetar model in order to match both the observed flux and time properties of the γ-ray signal. However, the optical light curve shows late-time bumpy structures that could either be the signature of multiple CSM shell interactions or evidence of central engine activity (perhaps magnetic reconfiguration, the Lense-Thirring effect, or other mechanisms). Plausible scenarios are therefore either a hybrid (magnetar+CSM) model or a pure magnetar model whereby both the γ-ray and optical properties are fully explained by a magnetar and an infalling accretion disk.

5. SN 2017egm stands out both photometrically and spectroscopically from the other events in our sample. It exhibits the steepest rise and reaches the highest peak luminosity. The presence of W-shaped O II absorption lines suggests a higher radiation temperature or excitation by nonthermal radiation from the inner magnetar wind nebula compared to the other SLSNe I in our sample, which lack these features. In contrast, the SLSNe II in the sample show narrow Balmer emission lines indicative of strong CSM interaction, suggesting a different powering mechanism. These distinctions may explain why γ-ray emission is only observed in SN 2017egm. In addition, it is worth noting that there can be a large diversity of ejecta masses in SLSNe, which can range from a few solar masses up to around 40 solar masses (Blanchard et al. 2020; Gomez et al. 2024), leading to large differences in the γ-ray opacities and time evolution. The γ-ray visibility is a competition between increasing transparency and decreasing pulsar spin-down energy, which implies that a low-ejecta-mass object would be favored for γ-ray detection as the needed transparency would be reached earlier.

6. As the tera-electronvolt γ-ray emission in the CSM model is strongly suppressed due to γ − γ absorption, a CTAO detection of a SLSN would strongly favor the magnetar scenario. In this regard, we simulated the possible detectability of a SN 2017egm-like object in the magnetar model and concluded that in 50 hours of CTAO-North observations, a maximum horizon of ∼140 Mpc for a 5σ detection is obtained, which limits the SLSN samples reachable by CTAO to the nearby objects. Given the low number of known nearby SLSNe, a precise rate is hard to estimate but it is likely on the order of a few per decade within ∼140 Mpc.

Acknowledgments

We thank the anonymous referee for the careful reading of the manuscript and for constructive and helpful criticism. This work made use of gammapy, numpy, scipy, astropy, and matplotlib libraries. FA wishes to thank Ismael Perez-Fournon, Matt Nicholl, and Morgan Fraser for enlightening discussions. The work on the Fermi-LAT analysis was supported by CNES. GMD thanks Weili Lin for generously providing further details on the CSM shells and model from Lin et al. (2023), and acknowledges support from the Beatriu de Pinós postdoctoral grant programme of the Departament de Recerca i Universitats de la Generalitat de Catalunya (2024 BP 00244). This work has been supported by the grant PID2024-155316NB-I00 funded by MICIU /AEI /10.13039/501100011033 / FEDER, UE and CSIC PIE 202350E189. This work was also supported by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M and also supported by MCIN with funding from European Union NextGeneration EU (PRTR-C17.I1). BDM acknowledges support from the National Science Foundation (grant AST-2406637). EC would like to thank the National Science Foundation (NSF) for their support made possible by the NSF grant AST-2506735. EC would also like to thank NASA and the Smithsonian Astrophysical Observatory (SAO) for their support via the Chandra X-ray Observatory (CXO) theory grant TM4-25003X. The Flatiron Institute is supported by the Simons Foundation. IV acknowledges support by the ETAg grant PRG2159. PJP is funded by the European Union (ERC, project number 101042299, TransPIre). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

References

Appendix A: Fermi-LAT extended data

The list of events with a likelihood probability of association with SN 2017egm corresponding to Fig. 4 is provided in Table A.1. Only events with probabilities larger than 50% according to the procedure described in Sect. 3.2 are included.

The Fermi-LAT data points of the SED shown in Fig. 5 (bottom panel) are presented in Table A.2. These were estimated with the analysis configuration described in Sect. 2.3 and in the Bayesian block time interval.

Appendix B: CSM model: Shock crossing times

The CSM model assumes a shock propagating through a series of shells with $V_{\mathrm{s}}=d{R}_{\mathrm{s}}/dt=V_{s,0}{\left[\frac{t}{1\mathrm{day}}\right]}^{m-1}$, where the expansion parameter is m = 0.75 and the initial shock velocity V_s, 0 = 9.8 × 10³ km/s (Lin et al. 2023). Table B.1 provides average densities as well as forward shock crossing times for each shell, derived from Extended Data Table 1 Lin et al. (2023). As in Sects. 3 and 5.2, T₀ is set to the discovery time (57896 MJD).

Appendix C: Estimate of γ − γ absorption in a CSM interaction scenario

For a first estimate of the impact that γ − γ absorption could have on the CSM model prediction for CTAO, we derived the expected τ_γγ following Sect. 2.3 in Fang et al. (2020). Assuming isotropy in the blackbody component, we estimated the photon density at the shock (n_opt) during the relevant evolution using the results from Table 3 in Zhu et al. (2023). The derived optical depth (τ_γγ ∼ n_optσ_γγR_s) for photons of 0.1, 1, and 10 TeV is represented in Fig. C.1. It is apparent that the medium would be optically thick at least until the end of the emission interval. To derive the absorption directly applicable to Fig. 8, we computed the weighted time average τ_γγ for the Bayesian block interval (Fig. C.1, bottom). As was expected, γ − γ absorption is not relevant for the LAT energy band but would prevent CTAO detection (i.e., at approximately tera-electronvolt energies). For the specific case of a SN 2017egm-like object, since no emission is expected for the CSM model after 250–300 days (see Fig. 7), a CTAO detection would support a magnetar-powered scenario.

Fig. C.1. Top: Temporal evolution of τγγ for photons of 0.1, 1, and 10 TeV. Vertical gray lines represent the Bayesian block interval identifying the LAT emission, while the horizontal dashed line (τγγ = 1) represents the transition to a transparent medium (τγγ ≲ 1). Shaded areas represent the error propagation from temperature and luminosity measurements in Zhu et al. (2023). Bottom: Absorption above 10 GeV, averaged over the Bayesian block time interval in Fig. C.1. The overall shape is a result of the peaked SN blackbody spectrum and the cross section σγγ.

All Tables

All Figures

	Fig. 1. Comparison of the r band absolute magnitude (left) and the pseudo-bolometric (using the gri bands; right) luminosity light curves for the objects of our samples. The light curves have been aligned to r band peak; see Sect. 2.2 for more details. Fewer time bins appear in the luminosity panel because of the unfulfilled requirement of simultaneous g − r − i observations in some bins.
In the text

Fig. 2. Luminosity light curves in the 100 MeV – 100 GeV energy range over 16 yrs for each SN of our sample from the Fermi launch to August 2024 with a time bin of 6 months. A flux point is shown when the TS > 4, otherwise upper limits at the 95% confidence level are reported. In all time bins, the spectral index of the tested source is fixed to 2. The SN discovery date is also indicated. For a comparison across the sample, the derived Fermi-LAT flux is transformed to luminosity using the distance indicated in Table 1.

In the text

Fig. 3. Left: comparison of the optical bolometric luminosity13 taken from Lin et al. (2023) and the γ-ray 100 MeV–100 GeV luminosity (this work) for SN 2017egm along with the TS in each 15-day time bin and with larger time blocks presented in Sect. 3.2 in black. Middle and right: γ-ray and optical luminosity observations and predictions from the magnetar model tailored to SN 2017egm by Vurm & Metzger (2021) for different nebula magnetization parameters (see Sect. 5.1). Hereafter the T0 value is taken at the discovery date 23 of March 2017.

In the text

Fig. 4. Photon angular distance to SN 2017egm optical position as a function of time. Using the best-fit model of the ROI, the probability that a given photon is associated with the SN is illustrated by colors, while the size indicates the photon energy. The γ-ray time interval defined by the Bayesian block algorithm is illustrated with the shaded area. Coverage limitations due to the failure of one of the solar array rotator drives on Fermi resulted in the loss of approximately 18 days of data in March and April 2018.

In the text

Fig. 5. Top: residual TS map in the 0.1–100 GeV energy range in the time interval defined by the Bayesian block algorithm (2017-07-05 to 2017-10-25). The plus symbols represent the sources from the 4FGL-DR4 catalog and the SN optical position. Bottom: γ-ray luminosity SED of SN 2017egm in the same time interval. The spectral predictions from the magnetar model are shown for different nebula magnetization parameter (see Sect. 5.1 for more details).

In the text

Fig. 6. Multiwavelength maps of the region around the giga-electronvolt excess: radio (VLA; left) and X-rays (XMM; right). The white marker represents the best-fit position of the giga-electronvolt excess, and the dotted curve is its localization error ellipse at a 95% confidence level. The cyan and green markers locate the nearby galaxies NGC 3191 and MCG+08-19-017, while the red and blue correspond to NVSS J101948+462332 and SDSS J101839.48+462315.8, respectively. White circles mark the quasars SDSS J101858.02+461812.3 and SDSS J101918.93+462301.4 discussed by Li et al. (2026).

In the text

Fig. 7. Left: energy conversion into cosmic rays as a function of time for different η values (Eq. (1)), and compared with the best fit obtained from Naima. Vertical lines represent the crossing times of the inner radius of each shell. Scenarios for ti = 1 d (green) and ti = tbr (black) are displayed, both assuming Ecutoff = 1 PeV. Right: γ-ray flux expected from the interaction of each shell (η = 0.05) compared with the Fermi-LAT flux (100 MeV–100 GeV). The arrow represents the Bayesian block interval of the emission peak, while the TS is shown as shaded green bars. Both panels display the relevant time interval for the CSM shell model.

In the text

Fig. 8. Top: Comparison of CSM and magnetar (magnetization εB = 10−6) spectral models of SN 2017egm at a distance of 135 Mpc together with Fermi-LAT (this work) and VERITAS (Acharyya et al. 2023) observations. In the CSM model, the predicted flux is shown with or without γ − γ absorption (see Appendix C), highlighting the fact that no emission is expected above 50 GeV. A simulation of 50 hours of CTAO-North observations between 50 GeV and 10 TeV is also shown for the magnetar model (see Sect. 7 for more details). For clarity only detection flux points (where TS > 4) are shown and upper limits are omitted. Using the CTAO simulation shown above, below we display the significance of a SN 2017egm-like source at different distances in the magnetar model only.

In the text

Fig. C.1. Top: Temporal evolution of τγγ for photons of 0.1, 1, and 10 TeV. Vertical gray lines represent the Bayesian block interval identifying the LAT emission, while the horizontal dashed line (τγγ = 1) represents the transition to a transparent medium (τγγ ≲ 1). Shaded areas represent the error propagation from temperature and luminosity measurements in Zhu et al. (2023). Bottom: Absorption above 10 GeV, averaged over the Bayesian block time interval in Fig. C.1. The overall shape is a result of the peaked SN blackbody spectrum and the cross section σγγ.

In the text