© The Authors 2026

1. Introduction

The search for astrophysical neutrino sources has been one of the most challenging pursuits in multimessenger astronomy during the past decade. Improvements to neutrino telescopes and multiwavelength programs carried out by the astronomical community have unveiled a plethora of interesting candidates. These sources span a range of objects extending from nearby galaxies (IceCube Collaboration 2022; Sommani et al. 2025) to peculiar and powerful active galaxies across cosmological distances (Kadler et al. 2016; IceCube Collaboration 2018; Garrappa et al. 2019; Franckowiak et al. 2020; Buson et al. 2023; Garrappa et al. 2024). Rare systems such as tidal disruption events (TDEs; Stein et al. 2021, Reusch et al. 2022, van Velzen et al. 2024), consisting of an interaction between a star and a massive black hole, have also been suggested as possible neutrino sources. The introduction of a real-time alert distribution program by the IceCube South Pole Neutrino Observatory in 2016 (Aartsen et al. 2017) significantly enhanced the number of apparent spatial and temporal coincidence cases of rare, high-energy (> 100 TeV) neutrinos with astrophysical sources observed at various wavelengths. However, since the detection of a diffuse high-energy neutrino flux in 2013 (IceCube Collaboration 2013) and a number of improvements in its characterization (Abbasi et al. 2024a), its origins remain a mystery.

Among the proposed contributors to the diffuse high-energy neutrino flux, supernova (SN) explosions within dense circumstellar material (CSM) represent a compelling source population. These include Type IIn SNe (SNe IIn), which are a subclass of SNe that exhibit strong, long-lasting emission lines, as well as a blue spectral continuum that most likely come from the interaction of the SN ejecta with CSM (Schlegel 1990; Smith 2017; Gal-Yam 2017). In such systems, shock regions are formed where fast ejecta encounter wind or shells around the progenitor star, driving diffusive shock acceleration of ions up to TeV–PeV energies. The dense CSM environment is an ideal target for those hadrons to produce pions via inelastic pp interactions, which, in turn, decay into high-energy neutrinos and gamma rays. After the shock breaks out into the optically thick wind (Ofek et al. 2010), a collisionless shock is expected to develop and multiTeV neutrinos are expected to be produced efficiently on timescales that range from hours to months (Katz et al. 2012, Murase et al. 2014, Murase 2018).

In addition, recent evidence has revealed that a large fraction of SNe IIn are located within compact distributions of CSM (Bruch et al. 2021, 2023). The resulting neutrino flux in these large number of sources could account for a significant fraction of the observed > 10 TeV neutrino background (Waxman et al. 2025). Because these dense environments are usually opaque to gamma-ray radiation, high-energy neutrinos are the best messenger to probe cosmic-ray acceleration at TeV-PeV energies. Pitik et al. (2022) investigated the possibility that the neutrino TDE candidate AT2019fdr could be a SN IIn instead. A stacking analysis using seven years of IceCube data and a sample of > 1000 core-collapse SNe (CCSNe; including types IIn, IIP, and stripped-envelope events) found no significant signal above the background (Abbasi et al. 2023). The derived 90% confidence limit, assuming a E^−2.5 spectrum, constrains the emission in high-energy neutrinos (10³–10⁵ GeV) from SNe IIn to ≲1.3 × 10⁴⁹ erg. This corresponds to a contribution of ≲34% to the diffuse neutrino flux, assuming the E^−2.5 spectral shape (Abbasi et al. 2023). Although these results do not identify a population of detected neutrino emitters, the limits obtained still leave room for the association of a few sources with single high-energy events.

On 2012 Marchs 30, a ∼2.7σ alert from the IceCube optical follow-up program (OFU; Franckowiak et al. 2009) consisting of two ∼TeV neutrinos detected within < 2 s and with a 1.32 deg angular separation was found to be coincident with the SN IIn PTF12csy (Aartsen et al. 2015). However, the a posteriori significance of this coincidence was estimated only at a 2.2σ level.

Lu et al. (2026) discussed the association of two SN IIn with events from the IceCube realtime stream of track-like high-energy neutrinos. They found a low (p ∼ 0.7%) probability that these two coincidences were observed by chance. We discuss one of these two coincidences, occurring between the neutrino event IceCube-250421A and SN2025cbj, in this work. Another example is the recent study on the association between an interacting SN Ibn and a high-energy neutrino (Stein et al. 2025), which has been estimated to have a probability of p ∼ 0.3% to occur by chance.

In this work, we discuss the association between the SN IIn SN2025cbj and the high-energy neutrino IceCube-250421A detected by IceCube. In Sect. 2, we present the detection of the high-energy neutrino event and how it was found in coincidence with SN2025cbj in our real-time follow-up program. In Sect. 3, we present the photometric and spectroscopic observations of SN2025cbj and explain how we derived useful physical quantities to characterize the SN. In Sect. 4, we discuss the multimessenger coincidence by means of resampling simulations. In Sect. 5, we estimate an expected neutrino rate for SN2025cbj based on the general model from Katz et al. (2012).

2. The high-energy neutrino event IceCube-250421A

On 2025 April 21, at 17:06:08.07 UTC, IceCube detected a track-like event with a 48.1% probability (i.e., signalness) of being astrophysical in origin (IceCube Collaboration 2025). The event belongs to the Bronze alert selection stream, whose average signalness is 30%. The best-fit reconstructed direction of the neutrino is at RA_J2000 = 240.91 (+2.62, −4.44) deg, Dec_J2000 = 28.67 (+1.70, −1.70) deg, with the errors referring to the extension of the 90% PSF containment. The relatively large errors in the reconstructed direction are due to the fact that the muon track was crossing a short amount of detector volume at one of the detector corners. The event has an estimated energy of 151.43 TeV, and its false alarm rate due to atmospheric background is 1.0254 events per year (IceCube Collaboration 2025).

As part of the Large Array Survey Telescope (LAST) follow-up program for multimessenger alerts and given the favorable observability of the neutrino arrival direction, we automatically triggered observations of the neutrino 90% error region at 22:19:54 UTC of the same day, 5.23 hours after the IceCube detection. The observations covered about 90% of the neutrino footprint reported in IceCube Collaboration (2025) by observing a single predefined LAST survey field. The LAST automatic transient-detection pipeline (Konno et al., in prep.) found one candidate transient optical counterpart in the observed neutrino footprint at RA_J2000 = 242.22 deg, Dec_J2000 = 27.0 deg (Garrappa et al. 2025a). At the time of the LAST observations, this transient was already known as SN2025cbj, a SN candidate classified as Type IIn at redshift z = 0.0675 (Sollerman et al. 2025). This transient was discovered by the Zwicky Transient Facility (ZTF, Bellm et al. 2019) on 2025 February 20 (about 60 days before the neutrino detection), located at ∼1.8 deg distance from the IceCube best-fit localization, within the 90% PSF containment.

Among known neutrino-source candidates, a number of gamma-ray blazars has been reported coincident with IceCube-250421A in Fermi-LAT realtime follow-up observations (Pfeiffer et al. 2025). They report the coincidence of five known gamma-ray sources. Among these, four are associated with blazars in the 4FGL-DR4 catalog (Ballet et al. 2023) and one is unassociated. In addition, a new gamma-ray source is reported for the first time and tentatively spatially associated to a known blazar detected in radio bands. However, none of these sources have been significantly detected in Fermi-LAT data integrated over timescales of days to months prior the neutrino detection. Despite the so-called gamma ray and neutrino connection in blazars is still not clear (e.g., Garrappa et al. 2024), this discourages the association of any of these sources with IceCube-250421A in the hypothesis unlike previous candidates (e.g., IceCube Collaboration 2018, Garrappa et al. 2019). In addition, we searched for TDE events in the Transient Name Server¹ (TNS) and the ZTF Bright Transients Survey (Fremling et al. 2020; Perley et al. 2020; Rehemtulla et al. 2024) and found no known TDE source coincident with the IceCube-250421A footprint.

3. Observations of SN2025cbj

3.1. Photometry

Here, we describe the photometric measurements of SN2025cbj. We also derive the peak and explosion time from the light curves available from various surveys.

3.1.1. LAST photometry

We observed the field of IceCube-250421A with LAST during its commissioning phase. The LAST system consists of an array of optical telescopes located at the Weizmann Astrophysical Observatory² in the Israeli Negev desert, designed to study the variable and transient sky (Ofek & Ben-Ami 2020, Ofek et al. 2023a, Ben-Ami et al. 2023). The first phase of the LAST system will consist of 72 Rowe-Ackermann Schmidt Astrograph f/2.2 telescopes with 28 cm diameter (40 already operative). Each telescope provides a field of view of 3.3 deg × 2.2 deg and a plate scale of 1.25 arcsec/pixel on camera. LAST observes in clear filter and the typical limiting magnitude (AB) of the system is 19.6 (21.0) in 20 s (20 × 20 s) exposures. The LAST data reduction pipeline and its products are described in Ofek et al. (2023b), while the image subtraction and transient detection pipeline is discussed in Konno et al. (in prep.).

On the night of April 21, 2025, the observations of the IceCube-250421A field were carried out under bright moon conditions, with an average limiting magnitude of 19.70 (AB) over a total of 19 epochs (Garrappa et al. 2025a). The average magnitude of SN2025cbj during the 19 epochs was m_AB = 18.53 ± 0.07 and no significant variability was observed. LAST monitored the source for nine additional nights. In addition, prior to the neutrino detection, LAST serendipitously observed the field containing SN2025cbj during sky-survey operations on April 4, 2025 (i.e., 17 days before the neutrino’s arrival). On this night, the source was detected at m_AB = 18.52 ± 0.06 (statistical error only), comparable with the brightness at the neutrino arrival. The total number of visits (20  ×  20 s exposures) observed by LAST for this object is 128. The average flux during each night is shown in the light curve of Fig. 1 and listed in Table 1. The photometric calibration of each image is obtained with the method described in Garrappa et al. (2025b).

Fig. 1. Upper: Photometric light curves from LAST (blue), ZTF-g (green), ZTF-r (red), DDOTI/OAN (orange) with temporal x-axis centered on the fitted peak time at MJD 60754.5. The magenta marker is the Swift-UVOT measurement in the UVW1 band. The dashed purple vertical line marks the arrival of IceCube-250421A. Lower: Pseudo-bolometric light curve in the (4700–6231 Å) range. The dot-dashed line shows the power-law fit of the early light curve decay after the peak (∝tα, α = −0.43).

3.1.2. ZTF photometry

We used photometric measurements from ZTF (Masci et al. 2019) retrieved from the Automatic Learning for the Rapid Classification of Events (ALeRCE) alert broker (Förster et al. 2021). The transient SN2025cbj (a.k.a. ZTF25aagbrpb) has been observed by ZTF with g and r filters during sky-survey operations. In Fig. 1, we show the ZTF light curves. We list the measurements in Table 2.

3.1.3. Other multiwavelength follow-up observations

After the LAST report of the possible multimessenger coincidence, the source was observed by DDOTI/OAN (Atri et al. 2025; Watson et al. 2016), which reported a flux level in the w (white) band of m_w = 18.53 ± 0.04, consistent with LAST measurements (included in Fig. 1).

The X-ray Telescope (XRT; Burrows et al. 2005) on board the Neil Gehrels Swift Observatory (Swift) observed SN2025cbj in a single target of opportunity (PI Veres) visit after the detection of IceCube-250421A. The observations were carried out on 2025 April 30 at 12:46 UTC for a total exposure of 2259 s. We used the Swift-XRT data products generator to reduce the data (Evans et al. 2020). No significant source was detected in Swift-XRT observations and we calculated a 3-σ upper limit in the 0.3–10 keV range of < 1.2 × 10⁻¹³ erg cm⁻² s⁻¹ (unabsorbed flux, using N_H = 4.1 × 10²⁰). This corresponds to an upper limit on the luminosity of L_XRT ≤ 1.4 × 10⁴² erg s⁻¹. The source was also observed during the same time with the Ultra-Violet/Optical Telescope (UVOT; Roming et al. 2005) on board Swift for a total of 2566.191 s exposure, using the UVW1 filter. We detected the source in UVOT data with apparent magnitude of m_W1 = 19.61 ± 0.09, corresponding to f(2634 Å) = (1.561 ± 0.129) × 10⁻² mJy.

As part of the realtime follow-up program of the Large Area Telescope (LAT; Atwood et al. 2009) on board the Fermi satellite, a search for gamma-ray emission at the localization of SN2025cbj was conducted in light of the multimessenger coincidence reported by LAST. Integrating observations of one month before the neutrino detection, no significant gamma-ray emission was detected (Pfeiffer et al. 2025). An upper limit on the flux was reported for a point source with power-law spectrum of index 2.0 at < 1.13 × 10⁻⁸ ph cm⁻² s⁻¹ (95% confidence limit). This corresponds to an upper limit on the gamma-ray luminosity in the 100 MeV – 1 TeV band of L_γ ≲ 1.8 × 10⁴⁴ erg s⁻¹. These additional multiwavelength observations are summarized in Table 3.

3.2. Optical spectroscopy and classification

SN2025cbj was classified as a SN IIn on March 5 2025 (Sollerman et al. 2025) based on a spectrum taken by the “SED” machine (SEDM; Ben-Ami et al. 2012; Blagorodnova et al. 2018) mounted on the P60 telescope in the Palomar observatory. The spectrum and classification were reported to the TNS and showed persistent narrow Balmer emission lines indicating long-lasting CSM interaction; thus it was classified as Type IIn.

Following the identification of the IceCube-250421A/SN2025cbj coincidence, we initiated a spectroscopic follow-up to confirm the SN was still interacting with CSM, which would allow a CSM-driven neutrino production channel. We initially obtained two spectra with the robotic SPRAT spectrograph (Piascik et al. 2014) mounted on the 2-meter Liverpool Telescope (LT) at the Roque de los Muchachos Observatory on La Palma (see Fig. 2). The spectra were obtained on April 25 and 28, 2025 and were reduced by the automatic SPRAT pipeline. We identified narrow Hα components in both spectra. However, due to the low brightness of SN2025cbj, we could not identify any other lines and, thus, we were not able to positively determine whether the narrow Hα emission was from the host galaxy, rather than from CSM interaction.

Fig. 2. Liverpool telescope spectra of SN2025cbj. A narrow Hα component can be identified in both epochs, below a strong galactic H α narrow feature. However, the lower resolution does not allow accurate velocity measurements or a positive identification of the Hβ line. We note that the narrow single pixel component in the 22.91 day spectrum is the result of instrumental noise.

Subsequently, we obtained a third, higher resolution (R ∼ 1340), spectrum using the BINOSPEC (Fabricant et al. 2019) spectrograph mounted on the 6.5-meter Multiple Mirror Telescope (MMT) at the MMT Observatory, Arizona on 2025 May 16. We present this spectrum in Fig. 3. The observation consisted of 3 × 450 s exposures. The data were acquired with a grating of 270 lines mm⁻¹ and a 1.0 ″ slit mask. The processed images, after the bias subtraction and flat-fielding using the Binospec pipeline (Kansky et al. 2019), are downloaded from the MMTO queue observation data archive. The spectra are reduced with IRAF, including cosmic-ray removal, wavelength calibration (using arc lamp frames taken immediately after the target observation), and relative flux calibration with an archived spectroscopic standard observation.

Fig. 3. MMT spectrum of SN2025cbj ∼56 days after the SN peak brightness. Zero velocity Hα and Hβ are marked with red and blue striped vertical lines, respectively. Zoom-ins on the Hα and β profiles with fits to a Lorentzian plus Gaussian shape are added, showing a good fit aside from the blue-wing asymmetry in Hα. The best-fit velocities (i.e., line FWHM) for both the Gaussian, vG, and Lorentzian, vL, components are presented next to the fits.

The MMT spectrum shows clear signs of CSM interaction, including narrow Hα and Hβ emission lines, superimposed on broader Lorentzian-shaped electron scattering wings, which are common in interacting SNe (Huang & Chevalier 2018). Interestingly, the Hα Lorentzian shape is very asymmetric, with its blue wing being much broader than its red wing. This is also seen to a lesser extent in Hβ. This effect has been observed in other interacting SNe such as SN 2013L (Andrews et al. 2017) and SN 2010jl (Smith et al. 2012; Ofek et al. 2014; Fransson et al. 2014) and indicates asymmetries in the CSM. However, since we only obtained a single epoch of medium-resolution spectrum, we cannot compare its evolution to that of other IIn SNe. Regardless, the consistency of narrow emission lines well after the neutrino detection confirms the SN had significant interaction at that time. All the spectra will be available on WISeREP (Yaron & Gal-Yam 2012) after the journal publication.

3.2.1. Narrow-line velocity

To measure the velocity of the narrow lines, we fit the Balmer lines in the MMT spectrum using the sum of a Lorentzian and a Gaussian profile, with the Lorentzian component representing the electron-scattering wings and the Gaussian the emission line. We then treated the full width half maximum (FWHM) of the Gaussian component as the bulk velocity of the emission line. We find both lines fit a velocity of v ∼ 230 km s⁻¹. However, this measurement should be treated as an upper limit, as the MMT spectral resolving power at the 270 lines/mm configuration is R ∼ 1340, which corresponds to a similar velocity limit of Δv ∼ 230 km s⁻¹. Regardless, the overall profile, including the much broader Lorentzian base, fits the lines well, indicating it is of CSM origin and not host galaxy emission lines. We also note that the Lorentzian line wings are too narrow to be of photospheric origin at v ∼ 3000 km s⁻¹, while typical SN photospheric velocities are closer to 10⁴ km s⁻¹ (Huang & Chevalier 2018), confirming the Lorentzian wings are indeed of CSM electron scattering origin. We present these fits as insets in Fig. 3.

3.3. Pseudo-bolometric light curve

We calculated a pseudo-bolometric optical light curve (4770-6231 Å) using observations in the ZTF-r and ZTF-g bands, which are the measurements that cover most of the light curve. Since ZTF observations are typically made with a single filter per night, we rebinned the light curve with a bin width of two days and considered only those bins containing at least one ZTF-r and one ZTF-g observation. The pseudo-bolometric luminosities were computed by fitting and then integrating a simple power-law model between the effective wavelengths of each band (λ_g = 4770 Å, λ_r = 6231 Å) and the pairs of fluxes in each two-day bin (F_g, F_r). The errors were derived via error propagation and assuming a 5% error on the measured distance. The pseudo-bolometric light curve is shown in the bottom panel of Fig. 1, the luminosity values are integrated only between the r and g bands, and the errors are statistical only. These values can be interpreted as lower limits on the bolometric luminosity.

3.4. Explosion and peak time

To estimate the explosion time (i.e., the time of first light), we used the earliest available photometric data for SN2025cbj, which is the ZTF-r data. The first r-band detection of m_g = 19.19 ± 0.1 was registered on 2025 February 19 at 12:07:02 UTC (MJD = 60725.5), following a nondetection on February 17, 2025 at 09:50:29 UTC (MJD = 60723.41) with a 5σ limiting magnitude of m_r > 19.5. We fit a broken power-law of f = a × (t − t_exp)ⁿ, where a is a fitting constant and t_exp is the explosion time, to the early r-band curve in flux space, while fixing the flux to 0 where t < t_exp. We fit all data up to the g-band maximum. We find the best fit to be a power law of n = 0.344 ± 0.056 and t_exp = 60723.21 ± 1.22 (best fit at February 17, 2025 at 05:14:17 UTC). While this result may seem contradictory to the original nondetection, preceding it by ∼4 hours, the error bar for our fit is quite large, and the nondetection is not very deep, indicating the explosion time could indeed precede it. We also note that SNe interacting with CSM have been shown to have a complex rise, which could change its trend very early, further increasing the uncertainty of the explosion time (as seen in the case of, e.g., SN 2023ixf; Li et al. 2024; Zimmerman et al. 2024).

We estimated the peak time of the ZTF light curve by fitting a simple Gaussian profile for both the ZTF-g and ZTF-r light curves. We find a best-fit value of t_peak, g = 60754.5 ± 1.4 MJD and t_peak, r = 60770.3 ± 3.7 MJD for the g and r bands, respectively. In the rest of this work, we use the value of t_peak, g = 60754.5 ± 1.4 MJD, which is better constrained by the observations.

3.5. Early light-curve decay

We fit the early decay of the pseudo-bolometric light curve after the peak time (MJD 60754.5) with a power-law function (∝t^α). We find the best-fit value at α = −0.43 ± 0.04. This value is consistent with expectations from interactions with massive CSM as well as with observations (e.g., SN2010jl, Svirski et al. 2012, Ofek et al. 2014). This supports the hypothesis of strong CSM interaction in SN2025cbj. The power-law fit of the early decay after the peak is shown in the bottom panel of Fig. 1.

We also fit the decaying part of the light curve with an exponential profile described by L = L₀ e^−t/τ_d, where t is the time since the peak and τ is the exponential timescale. We find the best-fit exponential timescale of the decay at τ_d = 142.43 ± 11.28 day. Although with a large error, this timescale is longer than the one expected from typical ⁵⁶Co decay (∼111 day). This longer decay timescale is similar to that observed in peculiar SNe such as 2010ij (Ofek et al. 2014; Fransson et al. 2014; Jencson et al. 2016). It supports the interpretation that the light curve of SN2025cbj is powered by interaction of the SN shock with the CSM.

4. Estimating the multimessenger coincidence

We investigated the possible association between SN2025cbj and IceCube-250421A by estimating the probability of chance coincidence between the two detections. To do so, we considered a collection of IceCube high-energy track-like events, following the one used in Garrappa et al. (2024). This collection includes archival and real-time events from the ICECAT-1 catalog (from 2011 to 2020, Abbasi et al. 2024b) and real-time events distributed through the GCN network from 2016 until the event IceCube-250708A. The initial neutrino sample consists of 391 track-like events. For the SN sample, we collected all transients from two catalogs: Transient Name Server (TNS, Gal-Yam 2021) and ZTF Bright Transients Survey (Fremling et al. 2020; Perley et al. 2020; Rehemtulla et al. 2024). From TNS, we selected 491 transients classified as a Type IIn SN discovered between 2013 May 9 and 2025 July 8. The ZTF-BTS catalog contains a smaller sample of 166 sources discovered in the same time period. All objects in the ZTF-BTS sample are also found in the TNS sample. We refer to the two statistical tests henceforth as T_tns (against the TNS sample) and T_bts (against the ZTF-BTS sample).

The purpose of the statistical tests described in this section is to estimate the chance coincidence probability of observing one or more neutrino-and-SN-IIn coincidences within δt = 61 days from the SN-IIn discovery date (i.e., the time delay observed between the SN2025cbj discovery and the detection of IceCube-250421A), since we would expect the neutrinos to be produced in the interaction between the SN ejecta and the CSM. To homogenize the SN and neutrino samples, additional selection cuts are required:

• Selecting neutrino events detected from 2013 May 09 (the discovery time of the first SN-IIn in the TNS sample) for T_tns;

• Selecting neutrino events detected from 2018 August 14 (the discovery time of the first SN-IIn in the ZTF-BTS sample) for T_bts;

• Selecting neutrinos and SN-IIn at high galactic latitudes (|b|> 20);

• Selecting neutrino events with a 90% PSF containment that is not larger than the one of IceCube-250421A (21.06 sq. deg).

The number of SNe and neutrino alerts for the final samples used in T_tns and T_bts are listed in Table 4. We simulated 10⁵ samples of neutrino alerts from the aforementioned list of IceCube neutrinos by scrambling their right ascensions while keeping their declinations and 90% error regions fixed. This preserves the specific distribution of IceCube alerts as a function of the arrival declination, which reflects the detector sensitivity across the sky. By repeating the search for coincident SNe IIn for each of the simulated samples within 61 days from their discovery date, we can empirically determine the spatial and temporal chance coincidence. We define a coincidence as an instance when there is at least one SN that is spatially coincident with the 90% error region of a neutrino alert, with the neutrino detection within a time δt from the SN discovery. Neutrinos coincident with multiple SNe are treated as a single coincidence.

4.1. Null hypothesis: the probability for random coincidences

From the statistical test T_tns, described in Sect. 4, we find that the distribution of spatial and temporal coincidences with SNe IIn in the background simulations follows a Poissonian distribution with best-fit λ = 0.273. Therefore, the probability of observing one or more coincidences for neutrinos with a reconstruction better than IceCube-250421A in T_tns is p(k ≥ 1) = 0.24. From the statistical test T_bts, we find that the number of spatial and temporal coincidences in the background simulations follows a Poissonian distribution with best-fit λ = 0.082, which gives p(k ≥ 1) = 0.078. The results are summarized in Table 4. These p-values are both consistent with the background hypothesis (i.e., the coincidence is observed only by chance). Regarding the two different SN IIn samples, the TNS database is a richer resource compared to the ZTF-BTS catalog because it includes contributions by a larger number of survey facilities.

4.2. Coincidences between SNe IIn and IceCube high-energy track events

In the previous section, we discuss the statistical significance of the SN-neutrino coincidence we are focusing on in this work. The purpose of this section is to have a more general view on similar coincidences that can be found in archival data. Since SNe are transient sources, it is key to define a temporal timescale in which to define the coincidence beside the spatial match. As we discuss in more detail in Sect. 5, in typical SNe-IIn interacting with CSM, a reasonable order of magnitude for the timescale of efficient neutrino production is of 𝒪(100) days. We ran the spatial and temporal search between the full sample of 343 IceCube alerts and 491 SN IIn from the TNS catalog. We list all the coincidences in Table 5. There are four high-energy neutrino alerts which are coincident with a SN IIn in TNS. One of these, IceCube-240307A, is coincident with two SNe IIn: SN2024bcm and SN2024cui, with δt = 39.8 and 20.5 days, respectively. At least one of these two coincidences is of random nature, given the very large error region of the neutrino localization (> 700 sq. deg). We noticed that all of these spatial coincidences are observed with poorly reconstructed neutrino events, with IceCube-250421A being the case with both the smallest error area (21.1 sq. deg) as well as the smallest angular separation between the SN and IceCube best-fit localization (δθ = 1.8 deg). We repeated the simulations described in Sect. 4 and we found that the probability of observing four or more random coincidences between the IceCube alerts and the TNS sample is p(k ≥ 4) = 0.25.

When we repeated the same matching between the IceCube alerts and the ZTF-BTS sample, we found only three coincidences, as the two objects SN2021jmm and SN2024bcm are not present in the ZTF-BTS catalog. The probability of observing three or more random coincidences between the IceCube alerts and the TNS sample is p(k ≥ 3) = 0.07. These p-values are consistent with the hypothesis that these coincidences are random in nature.

As a final control test, we used samples of SNe Ia from the TNS and the BTS catalogs that are not expected to be high-energy neutrino emitters (Murase 2018). Repeating T_tns and T_bts with δt = 100 days, we found 34 (15) neutrino alerts coincident with SNe Ia, over a sample of 11204 (6650) SNe. We estimated a chance coincidence of p_Ia, Ttns (k ≥ 34) = 0.46 and p_Ia, Tbts (k ≥ 15) = 0.55, which confirms the robustness of our statistical test.

4.3. Comparison with previous works

A recent work by Lu et al. (2026) investigates the coincidence between SNe IIn and IceCube high-energy neutrinos. In a broad archival cross-search between real-time high-energy track events from IceCube and the sample of classified SNe IIn from the ZTF-BTS catalog, they searched for spatial and temporal coincidences within 100 days from the SN discovery time and the arrival of neutrinos with a 90% containment region smaller than 30 sq. deg. The sizes of their final neutrino and SN-IIn samples are 138 and 163, respectively. They found, in addition to the IceCube-250421A and SN2025cbj coincidence, another coincidence between IceCube-231027A and SN2023syz. They reported a probability of finding two or more coincidences in the sample of p(k ≥ 2) = 0.0067.

We repeated our statistical test with the same selection cuts of Lu et al. (2026) and the ZTF-BTS sample. Based on these selection criteria, we found the additional coincidence between SN2023syz and IceCube-231027A that was excluded in our original search of events reconstructed better than IceCube-250421A. We found a probability of finding two or more random coincidences to be p(k ≥ 2) = 0.006, consistent with the one found in Lu et al. (2026). The p-value for k≥ 2 is intriguing. However, we remark that this is very sensitive to the preselection criteria. Moreover, neutrinos with large uncertainties in their reconstructions are more likely to change their containment regions and arrival directions in updated releases, while previously identified multimessenger coincidences might disappear (e.g., Zegarelli et al. 2025).

We note that when we repeated the cross-match with the same conditions of Lu et al. (2026), but with the TNS catalog, we found the same two coincidences (see Table 5, where events IC210503A and IC240307A end up removed by the cut of 30 sq. deg on the neutrino 90% area). However, after running the randomizations, we estimated a chance coincidence of observing a probability of finding two or more random coincidences to be p(k ≥ 2) = 0.067. This is about ten times higher than the previous p-value and it is likely due to the larger sample of SN-IIn in the TNS in the same time frame investigated by Lu et al. (2026); namely, with 381 instead of 169. The p-value is of the same order of magnitude of that computed in the test described in Sect. 4.2. We conclude that the fact that the p-value is so sensitive to the sample sizes weakens the neutrino-and-SN IIn connection. Moreover, it is sensitive to the cuts on the neutrino sample (e.g., error region size). Furthermore, the error regions of such events (especially those with large error contours) could be sensitive to the details of the reconstruction technique (e.g., Zegarelli et al. 2025). Therefore, such error regions are likely to be uncertain and it is difficult to quantify this uncertainty in the analysis. A larger number of neutrino events are also required to establish such a connection (see also Waxman et al. 2025).

5. Expected neutrino rate from SN2025cbj

Although we have shown that the IceCube-250421A/SN2025cbj coincidence is not at a statistically significant level (which is severely affected by the poor reconstruction of the neutrino event itself), we cannot simultaneously rule out the multimessenger association. Therefore, we can draw some simple considerations in the context of SN2025cbj being a neutrino source in light of the evidence of CSM interaction from our spectroscopic observations (see Sect. 3.2).

Unlike choked-jet scenarios (Mészáros & Waxman 2001; Murase & Ioka 2013; He et al. 2018; Zegarelli et al. 2024) in which the neutrino production is efficient only in an optically thick environment at the early stages of the SN explosion, the proton acceleration and interaction with the CSM are efficient after the so-called shock breakout on timescales that can last up to several weeks (see e.g., Katz et al. 2012). Assuming the dense CSM has a wind density profile (and noting that its radial extent is much larger than the breakout radius³), the profile can be simply parameterized by only the breakout radius, R_bo, (defined at optical depth of c/v_bo) and the shock velocity at that radius, v_bo,

$$\begin{array}{c}\hfill \rho (r)=\frac{c}{v_{\mathit{bo}}}\frac{m_p}{{\sigma }_T{R}_{\mathit{bo}}}{(r/R_{\mathit{bo}})}^{-2},\end{array}$$(1)

where c is the speed of light, m_p the proton mass, and σ_T the Thomson cross section (assuming the CSM is dominated by hydrogen). The rise time of the bolometric luminosity (Wasserman et al. 2025) can be estimated as

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{peak}}\simeq 1.2\frac{R_{\mathit{bo}}}{v_{\mathit{bo}}}log(c/v_{\mathit{bo}}).\end{array}$$(2)

This assumes that the shock is driven in the CSM by an initially polytropic stellar envelope that expands into the CSM. In this case, the shock decelerates slowly as v ∝ r^−0.1 (with only a weak dependence on the exact initial stellar density profile; e.g., Wasserman et al. 2025). This weak deceleration proceeds as long as the swept CSM mass is lower than the ejecta mass, which is consistent with the obtained result of M_bo ≈ 0.2 M_⊙, see below.

We can use Eq. (2) to constrain R_bo. We estimated the explosion time in Sect. 3.1.2 to be t_exp = 60723.21 ± 1.22 (MJD) and the peak time to be t_peak = 60754.5 ± 1.4 (MJD). This gives Δt_peak  ∼  31 days. The breakout velocity, v_bo, could not be measured directly from the spectrum. Photospheric velocities are typically measured by the blueshift absorption component of the broad H alpha P-cygni. However, there is no clear such absorption seen here (i.e., the feature to the blue of the wing is too narrow; therefore, it is likely noise). Therefore, we could not directly measure the shock velocity and assume a typical value of v_bo = 10⁹ cm/s (which is not expected to differ significantly). Plugging Δt_peak and v_bo in Eq. (2), we get R_bo ≃ 6.6 × 10¹⁴ cm. We integrated the CSM density profile (Eq. (1)) up to R_bo and we find the CSM mass contained up to R_bo is M_bo ≈ 0.2 M_⊙. Therefore, we do not expect the shock to decelerate significantly during propagation of a few breakout radii (as the swept CSM mass is much lower than the ejecta mass, which is at the level of least a few M_⊙ for typical SNe IIn).

Since the coverage of the light curve evolution is limited to the optical band, with only one observation in the NUV W1 band ∼40 days after the peak, we cannot accurately determine L_bol, p from observations; however, the estimated value L_bol, p = 2.3 × 10⁴³ is compatible with the expectation from a CSM breakout scenario, $L_{bol,p}\simeq {10}^{43}(\frac{R_{\mathit{bo}}}{{10}^{14}\mathrm{cm}})(\frac{v_{\mathit{bo}}}{{10}^9\mathrm{cm}/\mathrm{s}}{)}^2\mathrm{erg}{\mathrm{s}}^{-1}$ (Wasserman et al. 2025). In addition, the optical bands behavior also seems consistent with the expectation from the CSM breakout scenario, and the single X-ray upper limit obtained at 40 days past peak luminosity is consistent with the X-ray luminosity being smaller than 10% of the bolometric luminosity at that epoch (Wasserman et al. 2025).

The transition between the radiation-mediated shock (RMS) to the collisionless shock (CLS) occurs around 0.3 R_bo (Wasserman et al. 2025), while the proton acceleration and pion production are efficient approximately until the shock reaches ∼10 R_bo (Katz et al. 2012), with approximately constant neutrino luminosity during this period. The high-energy neutrino IceCube-250421A was detected ∼63 days after the explosion (which corresponds to ∼7 R_bo, including the weak shock deceleration), whereas the pion production should be efficient until ∼96 days from explosion (∼10 R_bo). The total nonthermal energy can be obtained following Katz et al. (2012):

$$\begin{array}{c}\hfill E_{\mathit{nt}}\sim {10}^{49}(\frac{R_{\mathit{bo}}}{{10}^{14}\mathrm{cm}}{)}^2(\frac{v_{\mathit{bo}}}{{10}^9\mathrm{cm}/\mathrm{s}}{)}^{-1}\mathrm{erg},\end{array}$$(3)

where we assume that the accelerated protons carry a 10% fraction of the post-shock energy. We estimated a total E_nt ∼ 4.3 × 10⁵⁰ erg. This is an order-of-magnitude estimation, which does not change significantly if the shock deceleration is included. Assuming an E² power-law spectrum for the protons with cutoffs at E_p, min = 10² GeV and E_p, max = 10⁸ GeV, we calculated the expected neutrino rate for a detector with the effective area of the IceCube Astrotrack Bronze stream⁴, at the declination band of SN2025cbj. Integrating over a time window of 96 days, we find N_νμ ∼ 1.5 × 10⁻³ expected muon neutrinos in the realtime stream event selection. This is consistent with the expected rates in other multimessenger associations with single high-energy neutrinos (e.g., Rodrigues et al. 2021, Omeliukh et al. 2025).

We repeated the same exercise for neutrinos starting from E_ν, min = 1 TeV, using the effective area used in the ten-year point source analysis, in the declination band of SN2025cbj (Aartsen et al. 2020; IceCube Collaboration 2021). We find an expected total number of N_νμ ∼ 1.7 × 10⁻² muon neutrinos within 96 days from the explosion. To estimate the detection rate, we need to multiply by the number of massive CSM SNe IIn, 𝒪(10²), observed during the IceCube operations, which gives 𝒪(0.1 − 1) expected detections. However, such detections will not be statistically significant. This estimate does not take into account the completeness of the catalogs used in the analysis. Finally, we briefly note that the nondetection in Fermi-LAT data and the corresponding upper limits are consistent with the discussion on the gamma-ray optical depth and luminosity in Katz et al. (2012) and Waxman et al. (2025).

6. Summary

We present the first broad set of optical photometric and spectroscopic observations of the SN IIn SN2025cbj, found in coincidence with the IceCube realtime high-energy neutrino event IceCube-250421A. The multimessenger association was first identified during the LAST neutrino follow-up program, although the SN was already discovered around 60 days before the neutrino detection. We present in this work for the first time high-resolution spectroscopic observations of SN2025cbj, which show clear evidence of CSM interaction at the time of the neutrino arrival. Although the optical spectrum indicates a promising environment for hadronic acceleration and neutrino production in a collisionless shock after the breakout, here we discuss how several factors challenge the significance of the multimessenger association. First, the poor reconstruction of the neutrino event weakens the spatial association. The ∼60-day delay affects the temporal coincidence, but it is consistent with the expected time frame in which neutrino production efficiency is high in these environments. We find that the spatial and temporal coincidence between IceCube-250421A and SN2025cbj is consistent with the background hypothesis of being a chance coincidence, with the smallest p-value of p(k ≥ 1) = 0.078 found against the ZTF-BTS catalog. We note that these p-values are sensitive to the size of the SNe and neutrino samples used for the statistical test.

We additionally show how a broader search for neutrino and SN-IIn coincidences in the entire sample of IceCube alerts reveals two additional high-energy track-like events coincident with SNe IIn within a 100-day window since the SN discovery. However, randomizations of these samples reveal that the background hypothesis cannot be rejected for this emerging population of multimessenger coincidences.

In conclusion, while we cannot rule out the SN2025cbj-and-IceCube-250421A association, we did compute the expected rate of neutrinos from SN2025cbj, assuming a typical scenario of proton acceleration and interaction in a collisionless shock. We found N_νμ ∼ 10⁻³ expected muon neutrinos within the time period of ∼96 days after the explosion in which SN2025cbj should, in principle, be an efficient neutrino source. This is consistent with the detection of a single high-energy neutrino event during this time. The expected rate is about a factor of 10 higher for neutrino energies as low as E_ν, min = 1 TeV. We encourage performing further dedicated point-source analysis of IceCube data from the direction of SN2025cbj to assess empirical constraints on these numbers. However, these estimates can vary significantly depending on the true characteristics of SN2025cbj, which could not be fully obtained from the available observations (e.g., see Pitik et al. 2022 for a detailed study on the SN+CSM model parameter space).

The multimessenger association between IceCube-250421A and SN2025cbj does not provide any independent evidence of high-energy neutrino production in SN IIn systems, but it can be considered among the population of possible neutrino sources identified by realtime multimessenger and multiwavelength observations. In the near future, larger numbers of detected neutrinos, a continuous improvement of the realtime alerts, and synergies among the follow-up programs across the whole electromagnetic spectrum will be key to detecting an emerging population of sources that are currently below the detection level. These searches will also benefit from the advent of future larger detectors (e.g., IceCube-Gen2, Aartsen et al. 2021) as well as new missions dedicated to the discovery and characterization of such transient events (e.g., ULTRASAT, Shvartzvald et al. 2024).

Data availability

All optical spectra from Sect. 3.2 are available on WISeREP⁵. Tables 1 and 2 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/708/A223.

Acknowledgments

We thank the anonymous referee for useful comments on the manuscript. We thank B. Katz, E. Waxman, and B. Zackay for useful discussions. We thank P. M. Veres for useful comments to the manuscript and for triggering the Swift ToO. S.G. is grateful for the support of the Koshland Family Foundation. E.O.O. is grateful for the support of grants from the Willner Family Leadership Institute, André Deloro Institute, Paul and Tina Gardner, The Norman E Alexander Family M Foundation ULTRASAT Data Center Fund, Israel Science Foundation, Israeli Ministry of Science, Minerva, BSF, BSF-transformative, NSF-BSF, Israel Council for Higher Education (VATAT), Sagol Weizmann-MIT, Yeda-Sela, and the Rosa and Emilio Segrè Research Award. This research is supported by the Israeli Council for Higher Education (CHE) via the Weizmann Data Science Research Center, and by a research grant from the Estate of Harry Schutzman. S.A.S. Is grateful for the Zuckerman Scholars fellowship. PC acknowledges the support from the Zhejiang provincial top-level research support program. Based on observations made with the Liverpool Telescope operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. CMC receives funding from UKRI grant numbers ST/X005933/1 and ST/W001934/1. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

References

All Tables

All Figures

Fig. 1. Upper: Photometric light curves from LAST (blue), ZTF-g (green), ZTF-r (red), DDOTI/OAN (orange) with temporal x-axis centered on the fitted peak time at MJD 60754.5. The magenta marker is the Swift-UVOT measurement in the UVW1 band. The dashed purple vertical line marks the arrival of IceCube-250421A. Lower: Pseudo-bolometric light curve in the (4700–6231 Å) range. The dot-dashed line shows the power-law fit of the early light curve decay after the peak (∝tα, α = −0.43).

In the text

	Fig. 2. Liverpool telescope spectra of SN2025cbj. A narrow Hα component can be identified in both epochs, below a strong galactic H α narrow feature. However, the lower resolution does not allow accurate velocity measurements or a positive identification of the Hβ line. We note that the narrow single pixel component in the 22.91 day spectrum is the result of instrumental noise.
In the text

Fig. 3. MMT spectrum of SN2025cbj ∼56 days after the SN peak brightness. Zero velocity Hα and Hβ are marked with red and blue striped vertical lines, respectively. Zoom-ins on the Hα and β profiles with fits to a Lorentzian plus Gaussian shape are added, showing a good fit aside from the blue-wing asymmetry in Hα. The best-fit velocities (i.e., line FWHM) for both the Gaussian, vG, and Lorentzian, vL, components are presented next to the fits.

In the text