Is the High-energy Neutrino Event IceCube-200530A Associated with a Hydrogen-rich Superluminous Supernova?

The Zwicky Transient Facility follow-up campaign of alerts released by the IceCube Neutrino Observatory has led to the likely identification of the transient AT2019fdr as the source of the neutrino event IC200530A. AT2019fdr was initially suggested to be a tidal disruption event in a Narrow-Line Seyfert 1 galaxy. However, the combination of its spectral properties, color evolution, and feature-rich light curve suggests that AT2019fdr may be a Type IIn superluminous supernova. In the latter scenario, IC200530A may have been produced via inelastic proton-proton collisions between the relativistic protons accelerated at the forward shock and the cold protons of the circumstellar medium. Here, we investigate this possibility and find that at most 4.6 × 10−2 muon neutrino and antineutrino events are expected to be detected by the IceCube Neutrino Observatory within 394 days of discovery in the case of excellent discrimination of the atmospheric background. After correcting for the Eddington bias, which occurs when a single cosmic neutrino event is adopted to infer the neutrino emission at the source, we conclude that IC200530A may originate from the hydrogen-rich superluminous supernova AT2019fdr.


INTRODUCTION
In 2013, the IceCube Collaboration reported the detection of a flux of high-energy neutrinos of astrophysical origin, marking the beginning of the high-energy neutrino astronomy era.Despite the growing number of high-energy neutrino events detected by the IceCube Neutrino Observatory, the sources of the cosmic neutrino flux remain to be unveiled (Aartsen et al. 2018a;Abbasi et al. 2021a,b;Aartsen et al. 2020).
The growing number of cosmic neutrino alerts has triggered follow-up searches for coincident detection of electro-magnetic radiation, see e.g.Abbasi et al. (2021c); Garrappa et al. (2019); Acciari et al. (2021).On October 1st 2019, the IceCube Collaboration reported the detection of a muon track neutrino of likely astrophysical origin, IC191001A.This event has been suggested to be the neutrino counterpart of the tidal distruption event (TDE) candidate AT2019dsg which was discovered by the Zwicky Transient Facility (ZTF) -see e.g.Stein et al. (2021); Bellm et al. (2019).Various theoretical models have been discussed to interpret this likely association (Winter & Lunardini 2021;Liu et al. 2020;Murase et al. 2020), however the jetted version of these models is being challenged by the most recent work on the radio properties of AT2019dsg (Cendes et al. 2021;Mohan et al. 2021;Matsumoto & Piran 2021;Matsumoto et al. 2021).
More recently, the follow-up campaign of IceCube neutrino alerts carried out by the ZTF Collaboration has led to another transient association.On May 31st 2020, IceCube Collaboration (2020); Reusch et al. (2020a) detected another muon track candidate (IC200530A), which was suggested to be associated with the optical transient AT2019fdr/ZTF19aatubsj 1 located at redshift z = 0.2666.The IC200530A event was detected ∼ 394 days after the discovery of the transient (hereafter identified with the onset of the shock breakout) and about 300 days after the peak of the electromagnetic emission.This neutrino event has a recon-structed neutrino energy of E ν 80 TeV and a signalness larger than 50% (IceCube Collaboration 2020; Reusch et al. 2020a;Stein 2021).
The intriguing coincidence of two IceCube neutrino events with two ZTF transient sources has triggered searches by the ANTARES Collaboration (Albert et al. 2021) and led to stringent upper limits on the neutrino emission from both sources.In addition, the Baikal-GVD Collaboration is currently investigating clusters of neutrino events detected along the same angular directions of both ZTF sources (Suvorova et al. 2021).
AT2019fdr is located close to the nucleus of its host galaxy and shows strong narrow hydrogen emission lines within its spectra.This led to the initial classification of AT2019fdr as either a flaring active galactic nucleus (AGN) in a Narrow-Line Seyfert 1 galaxy (Frederick et al. 2021), or a tidal disruption event (Chornock et al. 2019).This has resulted in interpretations of IC200530A being associated with an accreting black hole transient event (Stein 2021).However, Yan et al. (2019) proposed that AT2019fdr is a hydrogen-rich superluminous supernova (otherwise named superluminous supernova of Type IIn, SLSN IIn).Hydrogen rich SLSNe exhibit strong narrow Balmer emission lines within their spectra, but are more luminous than standard type IIn supernovae (SNe IIn), achieving luminosities typically with M −20 at peak brightness (Gal-Yam 2012;Smith 2014;Gal-Yam 2019).The narrow emission lines within SNe IIn are indicative of interaction between the SN ejecta with a dense shell of surrounding circumstellar material (CSM) in which kinetic energy is efficiently converted into thermal energy.The high luminosity of SLSNe IIn is thought to be the result of either an highly energetic explosion [with typical energies E ej O(10 51 -10 52 ) ergs], interaction with an unusually massive CSM (Moriya et al. 2018), or some combination of the two scenarios.
Proton acceleration, even beyond PeV energies, could take place in the proximity of the SLSN shock expanding in the dense CSM.The interaction of these protons with those of the shocked CSM may lead to copious neutrino emission (Murase et al. 2011;Katz et al. 2011;Murase et al. 2014;Cardillo et al. 2015;Zirakashvili & Ptuskin 2016;Petropoulou et al. 2016Petropoulou et al. , 2017;;Murase 2018).In this work, we investigate the possibility that IC200530A originates from AT2019fdr, under the framework that this transient is a SLSN IIn.
This paper is organized as follows.After introducing the main features characterizing AT2019fdr in Sec. 2, we outline the setup adopted to predict the neutrino signal in Sec. 3. Our findings are presented in Sec. 4 together with a discussion on the dependence of the neutrino signal on the uncertain parameters characteristic of AT2019fdr.A discussion on our findings and caveats is reported in Sec. 5, followed by our conclusions in Sec. 6.In addition, Appendix A discusses how the AT2019fdr parameter space is constrained by the observational constraints on AT2019fdr that we apply from neutrino and electromagnetic data.We investigate the temporal evolution of the maximum proton energy as a function of the model parameters considered for AT2019fdr in Appendix B.

AT2019FDR: A TYPE IIN SUPERLUMINOUS SUPERNOVA
AT2019fdr exhibits many properties compatible with those of other documented SLSNe IIn from the literature.Spectroscopically, the event shows intermediate-width (∼1000 km s −1 ) Balmer emission lines combined with narrow galaxy emission lines from the host, superimposed upon a blue continuum2 .The intermediate width Balmer emission features are characteristic of interacting core-collapse SNe (SNe IIn and SLSNe IIn), see e.g.Gal-Yam (2017); Moriya et al. (2018).Although these lines are also observed within nuclear transients (AGN flares and TDEs), the lack of intermediate components to the other host galaxy emission features (e.g.O III) disfavors the interpretation of this transient as an AGN flare.It is unlikely that these features mark AT2019fdr as a TDE, as these events generally exhibit much broader emission profiles than seen in AT2019fdr (typically line widths ∼ 10 4 km s −1 , Hung et al. 2017;Charalampopoulos et al. 2021).
The photometric behavior of AT2019fdr shows several features within the multi-band light curve, as displayed in Fig. 1, consistent with interaction-powered SNe.Although the slow rise time (∼80 days in the rest frame) and lengthy decline of the transient can be interpreted under each of the three potential paradigms suggested for its origin, the photometric evolution of AT2019fdr is not smooth.The light curve has a clear bump close to the peak (around 60 days from first light in the rest frame) alongside the beginning of an apparent re-brightening feature around 70 days after the optical peak.Episodes of re-brightening have been observed within some SNe IIn (e.g.Stritzinger et al. 2012;Nyholm et al. 2017) and are attributed to changes in the CSM density and variable progenitor mass-loss rates.
The late-time evolution of the transient (> 160 days from peak brightness) exhibits a slower decline than either Co 56decay (from a standard Ni 56 powered light curve) or the α = −5/3 power-law decline predicted by models of fallback accretion in TDEs (e.g.Rees 1988), but consistent with the range of typically slow declines in interacting SNe (Nyholm et al. 2020).AT2019fdr also exhibits a gradual reddening in color from peak to late times (transitioning from g -r ∼ −0.12-0.2mag), a property not predicted in nuclear transients, which shows an almost constant optical color in the majority of their light curves, but accordant with observations of normal SN IIn (Taddia et al. 2013).Finally, preexplosion variability is also not observed within the ZTF and ATLAS imaging (Yan et al. 2019), which disfavors an AGN origin.
Given the redshift of AT2019fdr, it is not possible to recover its complete rise in the ZTF photometry.However,  (Patterson et al. 2019), ATLAS (Tonry et al. 2018;Smith et al. 2020) and Swift (Gehrels et al. 2004).The detection epoch of IC200530A is marked as the black dashed vertical line and was observed ∼ 394 days after the first optical detection of the SN in the observer frame.We display the time from estimated shock breakout (t bo ), along the x-axis.
non-detections in the ZTF g-band prior to first light place the breakout epoch 6 days (5 rest frame days) before the first ZTF detection (see Fig. 1).Fitting a low order polynomial to the rise of the ZTF curve suggests that the start of the optical light curve coincides with these non-detections.Based on this, we assume the onset of the shock breakout at the first detection of MJD = 58606 ± 6 days.
We also note that AT2019fdr was not the only source suggested to be associated with the neutrino event IC200530A.AT2020lam 3 and AT2020lls 4 were also suggested to be possibly associated, as they were found within a 90.0%localization of the neutrino event (Reusch et al. 2020a).AT2020lam was classified using the Nordic Optical Telescope as a Type II SN located at z = 0.033 (Reusch et al. 2020c).However, the spectrum and light curve showed no evidence of CSM interaction, necessary for neutrino producing, leading Reusch et al. (2020c) to suggest that it was not associated with the neutrino event IC200530A.
AT2020lls was also classified using the Nordic Optical Telescope, but as a Type Ic SN located at z = 0.04106 that occurred ∼ 8 days prior to the detection of IC200530A (Reusch et al. 2020b).As this source did not show broad absorption features consistent with a subclass of Type Ic SN called Type Ic-BL, which are commonly associated with off-axis gammaray bursts or choked jets, Reusch et al. (2020b) suggested this was not associated with the neutrino event IC200530A.

MODEL SETUP
3 https://wis-tns.weizmann.ac.il/object/2020lam 4 https://wis-tns.weizmann.ac.il/object/2020lls In this section, we introduce the method adopted to compute the neutrino spectral energy distribution from AT2019fdr and its temporal evolution, as well as the properties of AT2019fdr useful to this purpose.Details on the estimation of the neutrino flux and event rate expected at Earth follow.

Spectral energy distributions of protons and neutrinos
We assume a spherical, steady and wind-like circumstellar medium (CSM) with solar composition ejected from the massive progenitor in the final stages of its evolution, as sketched in Fig. 2. We define its number density profile as where Ṁ is the stellar mass loss rate, v w the wind velocity, m = µm H , with µ = 1.3 being the mean molecular weight for a neutral gas of solar abundance, and R the distance to the stellar core.
The interaction of the stellar ejecta with the CSM leads to the formation of a forward shock (propagating in the CSM) and a reverse shock (propagating back into the stellar ejecta).Both the forward and reverse shocks could, in principle, contribute to the neutrino emission.Working under the assumption that the ejecta density profile decreases steeply (Chevalier & Fransson 2003), we neglect the contribution of the reverse shock since the forward shock is expected to predominantly contribute to the total energy dissipation rate and dominate the particle acceleration observed in SN remnants (e.g., Ellison et al. 2007;Patnaude & Fesen 2009;Schure et al. 2010;Slane et al. 2015;Sato et al. 2018;Suzuki et al. 2020).Hence, we focus on the neutrino emission from the forward shock for the sake of simplicity.
Following Chevalier (1982); Moriya et al. (2013), we assume that spherically symmetric SN ejecta of mass M ej and kinetic energy E k expand in the surrounding CSM homologously.The CSM extends up to an external radius R CSM (see Fig. 5).The outer ejecta density profile, which is relevant for the interactions leading to neutrino production, scales as n ej ∝ R −s , where we assume s = 10.The shocked SN ejecta and CSM form a thin dense shell because of efficient radiative cooling.Being the thickness of the thin shocked shell much smaller than its radius, one can describe its evolution through the radius R sh (t).In the ejecta dominated phase, namely in the phase in which most part of the ejecta is still freely expanding (i.e., when the mass of the ejecta is larger than the swept-up CSM mass), the shock radius is given by (Moriya et al. 2013;Chevalier & Fransson 2016): (2) with the corresponding shock velocity v sh = dR sh /dt.
Because of the high CSM density, the forward shock is initially expanding in a radiation dominated region, and particle acceleration is not efficient (Weaver 1976 The indigo line represents the forward shock that propagates radially outwards.The black dotted line marks the location of the deceleration radius of the ejecta (R dec ).The latter is located at radii smaller than R CSM (as in this sketch) for a relatively large CSM mass compared to the ejecta mass or larger than R CSM for very massive ejecta and rarefied CSM; see Eq. 4. For extremely large M CSM /M ej , it is possible that R dec < R bo .
Bromberg 2008; Katz et al. 2011;Murase et al. 2011).Efficient particle acceleration takes place at radii larger than that of the shock breakout (R bo ), where initially trapped photons are free to diffuse out to the photosphere; the shock breakout radius is computed by solving the following equation: where κ es ∼ 0.34 cm 2 g −1 (Pan et al. 2013) is the electron scattering opacity at solar abundances, and c is the speed of light.When the SN ejecta mass M ej becomes comparable to the swept-up mass from the CSM, the ejecta enters the CSMdominated phase.This transition happens at the deceleration radius Note that R dec may be located at radii smaller than R CSM as shown in Fig. 2, or larger than R CSM according to the relative ratio between M ej and M CSM (i.e., if M CSM > M ej , then R dec < R CSM and viceversa).Furthermore, for M CSM extremely large with respect to M ej , R dec can even be smaller than R bo .For R > R dec , the forward shock radius evolves as (Suzuki et al. 2020) . (5) where we have assumed adiabatic dynamical evolution for the sake of simplicity.At radii larger than R bo , diffusive shock acceleration of the incoming CSM protons takes place.Following Finke & Dermer (2012); Petropoulou et al. (2016), the proton injection rate for a wind density profile is where the parameter α dictates the radial dependence of the shock velocity (v sh ∝ R α ), it is α = −1/7 in the free expansion phase (R < R dec ) and α = −1/2 in the decelerating phase (R > R dec ).The fraction of the shocked thermal energy stored in relativistic protons is ε p , while H(x) = 1 for x > 0 and zero otherwise.We set the proton spectral index k = 2 and the minimum Lorentz factor of the accelerated protons γ p, min = 1.The maximum Lorentz factor of protons (γ p, max ) is obtained by requiring that the acceleration timescale in the Bohm limit, t acc ∼ 20γ p m p c 3 /3eBv 2 sh (Protheroe & Clay 2004), is shorter than the total cooling timescale for protons: t acc ≤ t p,cool .B = 32πε B m p v 2 sh n CSM is the magnetic field in the post-shock region, whose energy density is a fraction ε B of the post-shock thermal energy density U th = (9/8)m p v 2 sh n CSM .The latter is obtained by considering the Rankine-Hugoniot jump conditions across a strong non-relativistic shock with compression ratio approximately equal to 4.
The most relevant energy loss mechanisms for protons are inelastic pp collisions and the cooling due to adiabatic expansion of the shocked shell, hence t −1 p,cool = t −1 pp + t −1 ad , with t pp = (4k pp σ pp n CSM c) −1 , where we assume constant inelasticity k pp = 0.5 and energy-dependent cross-section σ pp (E p ) (Zyla et al. 2020).Following Fang et al. (2020), the adiabatic cooling is t ad = min[t dyn , t cool ], where t cool is the typical cooling time of the thermal gas behind the shock and t dyn is the dynamical time of the shock.When the shock is radiative, the particle acceleration region is shrank to a characteristic length ∼ v sh t cool , limiting the maximum achievable particle energy.The cooling time is t cool = 3k B T/2n sh Λ(T ) (Franco et al. 1992) where k B is the Boltzmann constant, n sh = 4n CSM is the density of the shocked region, and Λ(T ) is the cooling function capturing the physics of radiative cooling.Here T is the gas temperature immediately behind the forward shock front obtained by the Rankine-Hugoniot conditions, given by: where γ = 5/3 is the adiabatic index of the gas.Finally, the cooling function [in units of erg cm 3 s −1 ] is given by the following approximation (Chevalier & Fransson 1994): 6.2 × 10 −19 T −0.6 105 < T 4.7 × 10 7 K 2.5 × 10 −27 T 0.2 T > 4.7 × 10 7 K .
(8) where line emission dominates at low T and free-free emission at high T .
Relativistic protons in the shocked region may also interact with the ambient photons via pγ interactions.However, in this work we ignore this energy loss channel, consistent with the work of Murase et al. (2011); Fang et al. (2020), which show that pγ interactions can be safely neglected for a wide range of parameters.
Since we aim to compute the neutrino emission, we track the temporal evolution of the proton distribution in the shocked region between the shock breakout radius R bo and the outer radius R CSM .
The evolution of the proton distribution is given by (Sturner et al. 1997;Finke & Dermer 2012;Petropoulou et al. 2016): where N p (γ p , R) represents the total number of protons in the shell at a given radius R with Lorentz factor between γ p and γ p + dγ p .The radius R is related to the time t measured by an observer at Earth: t = t(R)(1 + z), where we denote with a tilde all parameters in the reference frame of the central compact object hereafter.The second term on the left hand side of Eq. 9 takes into account energy losses due to the adiabatic expansion of the SN shell, while pp collisions are treated as an escape term (Sturner et al. 1997).Other energy loss channels for protons are negligible (Murase et al. 2011).Furthermore, in Eq. 9 the diffusion term has been neglected since the shell is assumed to be homogeneous.
The neutrino production rates, for muon and electron flavor (anti)neutrinos are given by (Kelner et al. 2006): ν µ , F (2) ν µ and F ν e follow the definitions in Kelner et al. (2006).Equations 10 and 11 are valid for E p > 0.1 TeV, corresponding to the energy range under investigation.

Parameters characteristic of AT 2019fdr
Numerical simulations aiming to model SLSNe IIn light curves are undergoing, see e.g.Dessart et al. (2015); Chatzopoulos & Tuminello (2019); Suzuki et al. (2021); Suzuki et al. (2019); however, the exact underlying physics which powers these sources is still uncertain.In the following, we outline the properties of AT2019fdr useful to model its neutrino emission.
By relying on existing data on AT2019fdr from ZTF (Patterson et al. 2019), ATLAS (Tonry et al. 2018;Smith et al. 2020) and Swift (Gehrels et al. 2004), we integrate the observed pseudo-bolometric light curve and estimate that the total radiated energy from AT2019fdr is Ẽrad = 1.66 ± 0.01 × 10 52 erg.To take into account the uncertainties on the radiative efficiency, namely the fraction of the total energy that is radiated, we consider two characteristic values of the kinetic energy Ẽk of the ejecta: 5 × 10 52 erg and 10 53 erg, which correspond to a radiative efficiency of ∼ 35% and 18%, respectively (see Chevalier & Irwin (2011), where the total radiated energy is estimated to be E rad = 0.32E k ).
We assume the proton fraction equal to ε p = 0.1 (Murase et al. 2011).This value is consistent with simulations of particle acceleration and magnetic field amplification at non-relativistic quasi-parallel shocks (Caprioli & Spitkovsky 2014).A discussion on the impact of different values of ε p on the expected neutrino event rate is reported in Sec. 5.The fraction of the post-shock internal energy that goes into amplification of the magnetic field is instead assumed to be The wind velocity is considered to be v w = 100 km s −1 (Moriya et al. 2014).The average mass loss rate is given by (Suzuki et al. 2021): where M CSM is the CSM mass contained within a shell of radius R CSM .
By fitting a basic polynomial to the bolometric light curve and available non-detections of AT2019fdr to extrapolate beyond the detection limits of ZTF, we estimate that its rise time (i.e., the time during which the luminosity reaches peak value, see Fig. 1) is t rise ∼ 98 days.In addition, in order to link t rise to the other model parameters characteristic of AT2019fdr, we rely on the following relation (Suzuki et al. 2020): where the diffusion time is the time required for the radiation to travel from R to R ph 5 , and τ T (R) is the optical depth of the CSM at radius R. The rise time is expected to increase as a function of M CSM , since a massive and dense CSM prolong the photon diffusion timescale.Yet, in order to predict the correct behavior of t rise , one should take into account the effect of the variation of all the parameters: E k , M ej , M CSM , and R CSM .
The exact values of M ej , M CSM , and R CSM are highly uncertain for AT2019fdr and degeneracies may be at play when interpreting the AT2019fdr light curve.The reprocessing of information from the explosion by interaction with the CSM masks the properties of the SLSN explosion underneath it.Although the CSM density can be estimated in several ways, e.g. from the strength of the H-α line (Taddia et al. 2013) or through X-ray and radio observations (Chandra 2018), AT2019fdr lacks the necessary time series multiwavelength and spectroscopic data required to constrain it.Hence, we consider ranges of variability for the most uncertain parameters: M ej ∈ [5, 150] M , M CSM ∈ [5, 150] M , and R CSM ∈ [2, 4] × 10 16 cm.Out of these, we only take into account the ones in agreement with the measured t rise (allowing for an uncertainty of 50%) and requiring that the production of the neutrinos observed by the IceCube Observatory at ∼ 394 days after the breakout takes place inside the CSM, namely t(R CSM ) − t(R bo ) 394 days.See Appendix A for more details.A summary of the default values for the parameters considered for AT2019fdr is reported in Table 1.The neutrino and antineutrino flux (F ν α +ν α with α = e, µ, τ) at Earth from a SN at redshift z and as a function of time in the observer frame is [GeV −1 s −1 cm −2 ]: where Q ν β +ν β is defined as in Eqs. 10 and 11.Neutrinos change their flavor while propagating, hence the flavor transition probabilities are given by (Anchordoqui et al. 2014): P ν e →ν µ = P ν µ →ν e = P ν e →ν τ = 1 4 sin 2 2θ 12 , (15) with θ 12 33.5 deg (Esteban et al. 2020), and P ν β →ν α = P νβ →ν α .The luminosity distance d L (z) is defined in a flat ΛCDM cosmology as where Ω M = 0.315, Ω Λ = 0.685 and the Hubble constant is The neutrino fluence [GeV −1 cm −2 ] is calculated using with t bo = t(R bo ) and the time integral being restricted to 394 days.Finally, the event rate of muon neutrinos and antineutrinos expected at the IceCube Neutrino Observatory is where A eff (E ν , δ) is the detector effective area (Abbasi et al. 2021a).The minimum neutrino energy is E ν, min = 100 GeV for the declination of interest (Abbasi et al. 2021a), and F ν µ +ν µ (E ν , t) has been introduced in Eq. 14.In the following, we work under the assumption of perfect discrimination between astrophysical and atmospheric neutrinos; see Sec. 5 for a discussion on the expected event rate if the event sample should be contaminated by atmospheric neutrinos in the energy region below 100 TeV.The maximum neutrino energy E ν,max is related to the maximum proton energy: E ν,max = xE p,max .The total number of muon neutrinos and antineutrinos is computed over the temporal interval of 394 days: 4. FORECAST OF THE NEUTRINO SIGNAL In this section, we present the results on the neutrino signal expected from AT2019fdr.First, we discuss the neutrino spectral energy distribution and the event rate expected in the IceCube Neutrino Observatory.We then investigate the dependence of the expected signal on the uncertainties of the SLSN IIn model.

Energy fluence and temporal evolution of the neutrino event rate
Before focuing on the energy fluence and event rate of the detectable neutrino signal, we explore the characteristic cooling times of protons and the acceleration timescale characteristic of AT2019fdr, introduced in Sec.3.1.In order to give an idea of the variation of the cooling and acceleration timescales across the SLSN shell, Fig. 3 shows the proton cooling times as a function of the proton energy in the reference frame of the central compact object and at the representative radii R bo and R CSM for the SLSN configuration with ( Ẽk , R CSM , M ej , M CSM ) = (10 53 erg, 4 × 10 16 cm, 6 M , 49 M ).As discussed in the following, this SLSN configuration leads to the most optimistic scenario for neutrino production.
Proton-proton collisions are responsible for the dominant energy loss channel.Even though Fig. 3 represents the characteristic cooling times for one specific SLSN configuration, the hierarchy between pp and adiabatic losses is representative of all SLSN configurations considered in this work (lower Ẽk and R CSM larger than the ones adopted here would lead to scenarios with adiabatic energy losses being dominant over pp ones).
The evolution of E p,max depends on the specific choice of parameters Ẽk , R CSM , M ej , and R CSM , determining whether R bo ≶ R dec .For the typical values of Ẽk and R CSM considered in this work, the condition t pp < t ad is always fulfilled, and E p,max increases as a function of R up to R dec , and decreases otherwise.In fact, by using Eqs. 1, 2 and 5, we find: (22) Appendix B provides more details on the scaling of E p,max as a function of the SLSN model parameters.
The muon neutrino and antineutrino fluence, defined as in Eq. 19, is shown in Fig. 4 as a function of the neutrino energy.The band takes into account the uncertainties on the parameters characterizing AT2019fdr (see Sec. 3.2) and is defined by the parameter configurations leading to the highest and lowest neutrino fluence.Note that, for the SLSN parameters adopted in this work, the synchrotron cooling of charged pions and muons produced via pp interactions is negligible.In fact, the typical energies for which this energy loss becomes relevant are at least three orders of magnitude larger than the maximum achievable proton energies.Therefore, the neutrino spectra are not affected by the cooling of mesons.
Given our selection criterion (i.e., the observation of IC200530 about 394 days after the shock breakout and the constraints on the rising time of the light curve of AT2019fdr), the scenarios with the lowest fluence are the ones corresponding to configurations with large R CSM , low M CSM and high M ej .On the other hand, given the reduced parameter space allowed for low R CSM (see Appendix A), the most optimistic scenario corresponds to the highest R CSM , the lowest accessible M ej , and intermediate values of M CSM (M CSM 30-50M ).We refer the reader to Sec. 4.2 for a discussion on the dependence of the neutrino fluence from the SLSN characteristic parameters.
The reconstructed neutrino energy for the IC200530 neutrino event is marked with a dotted vertical line and it falls in the same energy range as the predicted fluence.One can see that, around the reconstructed energy of IC200530, the fluence can vary up to O(10 5 ) in magnitude.However, it is worth noting that the reconstructed energy carries an intrinsic uncertainty and may differ from the real energy of the detected neutrino, nevertheless we show it in order to guide the eye.
The muon neutrino and antineutrino event rate expected in IceCube (Eq.20) is shown in Fig. 5 as a function of time.The band in Fig. 5 takes into account the uncertainties on the characteristic quantities of AT2019fdr summarized in Table 1.For all SLSN cases within the envelope in Fig. 5, the event rate increases rapidly at early times.After the peak, depending on whether R dec > R bo or R dec < R bo , the neutrino event rate has a steeper or shallower decay.These two different trends are related to the evolution of the shock velocity and the maximum proton energy E p,max .Indeed, E p,max increases up to R dec as t increases and declines later.Since the detector effective area A eff increases as a function of E ν (Abbasi et al. 2021a) and the decline of v sh for R bo < R < R dec is shallow, a compensation effect can arise among the two quantities; hence, the drop of the Ṅν µ +ν µ curve can be slow.Viceversa, when both E p,max and v sh decrease, the event rate drops faster.Around the day of detection of IC200530 (marked by It is important to note that only a sub-sample of the SLSN parameter set reported in Table 1 allows us to obtain a neutrino signal compatible with our observational constraints.For example, none of the SLSN scenarios with Ẽk = 10 53 erg and R CSM = 2 × 10 16 cm passes our selection criteria, since the shock crosses the CSM envelope in a time shorter than 394 days.

Dependence of the neutrino signal on the parameters of AT2019fdr
In order to better explore the dependence of the neutrino signal expected in IceCube on M ej and M CSM , for Ẽk = 10 53 erg, first we investigate the neutrino fluence as a function of M CSM for fixed R CSM and M ej and then we fix M CSM and vary M ej .The choice of M CSM and M ej is guided by the SLSN configurations that better highlight the changes in the neutrino fluence for R bo ≶ R dec .From the panel on the left in Fig. 6, we see that the fluence increases as M CSM increases up to M CSM = 85 M .For larger M CSM , R bo > R dec , and therefore a turnover with a slow drop can be observed.Furthermore, a slight shift of the neutrino cutoff energy towards lower energies is visible as M CSM increases.The latter is due to the enhanced pp energy loss determined by the larger den-  1.The event rate increases rapidly at early times.After peak, the event rates for the SLSN scenarios representing the edges of the envelope decline because of the dominant decreasing trend of v sh as a function of time.In some intermediate scenarios, the increasing trend of E p,max and shallow decrease of v sh can be compensated, providing an increasing event rate at the moment of the detection.The neutrino event IC200530 has been observed ∼ 394 days after t bo as indicated by the dotted vertical line.In the proximity of the detection day, the event rate can vary up to a factor O(10 3 ) in magnitude.sity as well as the smaller v sh , which prevent particles from being accelerated to higher energies (see Eq.22).
In the right panel of Fig. 6, we observe an enhancement of the fluence as M ej decreases.Nevertheless, this trend is inverted for M ej 13 M , representative of the regime with R bo > R dec , where the lower v sh is responsible for a slight decrease in the neutrino production, together with a shift of the neutrino energy cutoff to lower energies.
Figure 7 shows the temporal evolution of the muon neutrino and antineutrino flux for the scenarios with the highest (left panel) and the lowest (right panel) expected number of neutrinos.In all cases, the flux decreases as time increases and shifts to lower or higher energies, for the most optimistic and pessimistic scenarios, respectively.Around the day of detection, the flux in the best scenario is a factor O(10 5 ) larger than the most pessimistic scenario.
In order to investigate the origin of IC200530, we integrate the event rate over 394 days of the neutrino signal for all selected SLSN configurations and obtain the total number of muon neutrino and antineutrino events, N ν µ +ν µ (Eq.21).A contour plot of N ν µ +ν µ in the plane spanned by M ej and M CSM is shown in Fig. 8 for R CSM = 4 × 10 16 cm and Ẽk = 10 53 erg as a representative example.The allowed region of the pa- The black solid lines marks the allowed region of the parameter space, defined by requiring that the location of the shock at the day of neutrino production is still in the CSM envelope and that the SLSN model parameters are compatible with the the light curve rise time.For fixed M CSM , the total neutrino number decreases as M ej increases, given the decline of the shock velocity v sh .Viceversa, for fixed M ej , the number increases as M CSM increases, given the larger number of proton targets for pp interactions.In the region R bo > R dec , one can see an inverted trend.The dotted lines correspond to the contour lines which track the scenarios providing the number of neutrino events displayed, and show how the neutrino number is affected in the transtition from R bo > R dec to R bo < R dec regions.See the main text for more details.

DISCUSSION
Table 2 summarizes the total number of muon neutrino and antineutrino events expected within 394 days from the shock breakout from AT2019fdr for the most optimistic and pessimistic SLSN configurations in terms of neutrino emission.The largest [smallest] number of events is obtained for the SLSN configuration with ( Ẽk , R CSM , M ej , M CSM ) = (10 53 erg, 4 × 10 16 cm, 6 M , 49 M ) [(5 × 10 52 erg, 4 × 10 16 cm, 150 M , 19 M )], and correspond to the edges of the band in Fig. 5.
An important aspect to consider in the interpretation of the neutrino event IC200530 concerns the discrimination of the atmospheric neutrino background, which dominates over the astrophysical neutrino flux below 100 TeV.As such, in Table 2 we distinguish between one case with the lower energy cutoff fixed at 100 GeV, mimicking excellent discrimination of the atmospheric neutrino background (see Sec. 3.3), and one more conservative case with the lower Table 2. Number of muon neutrino and antineutrino events expected in 394 days from the shock breakout from AT2019fdr for the most optimistic and pessimistic scenarios, with the low energy cutoff fixed at 100 GeV (i.e., excellent discrimination between the astrophysical and atmospheric signals) and 100 TeV (i.e., under the conservative assumption that the atmospheric background could not be eliminated).The most optimistic and pessimistic scenarios correspond to the following SLSN model parameters: ( Ẽk , R CSM , M ej , M CSM ) = (10 53 erg, 4 × 10 16 cm, 6 M , 49 M ) and (5 × 10 52 erg, 4 × 10 16 cm, 150 M , 19 M ), respectively.In the last column we estimate the signalness [N νµ+νµ,astro /(N νµ+νµ,astro + N νµ+νµ,atm )], by computing the number of atmospheric neutrino events over a period of 394 days, for 0.75 • around the direction of the source.energy cutoff at 100 TeV.The latter case reproduces a situation where the atmospheric neutrino events could not be distinguished from the astrophysical ones in the lower energy range.Although a dedicate likelyhood analysis is beyond the scope of this work, the last column of Table 2 reports N ν µ +ν µ ,astro /(N ν µ +ν µ ,astro + N ν µ +ν µ ,atm ), which should give an idea of the expected signalness and gives an indication of the probability that a detected neutrino event could be of astrophysical origin.It is evident that only an optimal discrimination of the atmospheric neutrino background allows to obtain a signalness of 40%, roughly comparable with the one of the neutrino event IC200530.The evolution of the neutrino curve shown in Fig. 5 should be considered carefully.In fact, some intermediate SLSN scenarios enclosed in the envelope in Fig. 5, and compatible with the reconstructed energy of the neutrino event IC200530A, have an event rate still increasing at the day of detection, therefore increasing the neutrino detection chances at later times, as it is the case for the neutrino event IC200530.
In order to assess whether the number of expected events in Table 2 is compatible with the detection of one neutrino event (IC200530) after 394 days from the shock breakout, we take into account the Eddington bias on neutrino observations.The Eddington bias must be taken into account when dealing with very small number of neutrino events, such as in this case; we refer the interest reader to Strotjohann et al. (2019) for a dedicated discussion.By relying on the local rate of SLSN IIn provided in Quimby et al. (2013) and integrating over the cosmic history by assuming that the redshift evolution of SLSN IIn follows the star formation rate (Yuksel et al. 2008), we obtain an average effective density of SLSN IIn equal to O(3 × 10 3 ) Mpc −3 .Although Fig. 2 of Strotjohann et al. (2019) was derived within a simplified framework and for constant redshift evolution, by extrapolating to larger effective source densities we conclude that the number of expected events in Table 2 may be compatible with the detection of at least one or two neutrino events from AT2019fdr.By taking into account the fact that the neutrino energy distribution of AT2019fdr falls in a region where the discrimination of the atmospheric neutrino background may be challenging, our findings hint towards a possible association of the neutrino event IC200530 to AT2019fdr.In addition, our results are compatible with the upper limits on the neutrino emission from the AT2019fdr source placed by the ANTARES Collaboration (Albert et al. 2021).
We should stress that the forecasted number of expected neutrino events includes some caveats related to our modeling.For example, one of the sources of uncertainty in the computation of the neutrino flux is the proton acceleration efficiency ε p .In this paper, we have adopted an optimistic ε p = 0.1, assuming that the shocks accelerating protons are parallel or quasi-parallel and therefore efficient diffusive shock acceleration occurs (Caprioli & Spitkovsky 2014).However, lower values of ε p would be possible for oblique shocks, with poorer particle acceleration efficiency.Values as low as ε p 0.003-0.01have been inferred in Aydi et al. (2020) for a nova, assuming shocks as the powering source of the simultaneously observed optical and γ-rays.However, observational constraints from other optical transients, including SLSNe, are still lacking; in addition, AT2019fdr is much more luminous than classical novae, possibly hinting to different conditions present in the acceleration region.
We stress that the neutrino flux scales linearly with ε p , allowing the reader to easily scale our results.All cases summarized in Table 2 would be compatible with the detection of one neutrino event, after taking into account the Eddington bias.Indeed, the detection of a single neutrino event may actually hint towards intermediate SLSN configurations, as well as values of ε p lower than our benchmark one.
Similarly, in this work we have assumed that protons are accelerated at the shock to a power law with slope k = 2, which is predicted by the test particle theory of diffusive shock acceleration.Nonetheless, non-linear effects involving the amplified magnetic field can kick in, modifying the shock structure and making the cosmic ray spectra mildly steeper than k = 2 (Caprioli et al. 2021).Larger k would result in steeper neutrino spectra, since the latter inherit the shape of the parent proton spectrum; as a consequence, lower fluxes should be expected in the energy of interest.
Another caveat to take into account concerns the use of the AT2019fdr discovery date in the observer frame as the breakout time of the shock.In fact, based on the non-detections in the ZTF data, we have assumed an explosion epoch at the first detection at MJD = 58606 ± 6 days on the basis of a fit on the existing data.Nevertheless, even allowing for an onset of the shock breakout to be as much as ∼ 20 days earlier than the first observed light, our predictions in Table 2 would not be affected by a factor larger than 10%.
Since initial submission of this manuscript, other publications have analysed IC200530 under the paradigm of a TDE origin (Reusch et al. 2021).The additional data presented within these works suggest that an apparent increase in the late time near infrared (NIR) emission may be attributed to a dust from the TDE emission.However, increasing late time NIR emission has been seen in other interacting SNe.For instance, the bright SN IIn SN2010jl exhibits such a NIR increase at late times; high-resolution spectroscopic observations show that this increasing emission was the result of rapid dust formation within the SN ejecta (Gall et al. 2014).
In addition, the vast majority of TDEs show bright X-ray emission over the full optical/UV evolution of an event (e.g., Auchettl et al. 2017;Brown et al. 2017;Hinkle et al. 2021;Wevers et al. 2021).Of those whose emission is dominated by optical/UV but has been detected in X-rays, the X-ray luminosities are an order of magnitude or more fainter than the eROSITA detection (e.g., Jonker et al. 2020;Holoien et al. 2019;Hung et al. 2020Hung et al. , 2021)).In addition, AT2019fdr is found close to the nucleus in a Narrow-Line Seyfert 1 active galaxy (Frederick et al. 2021).Seyfert AGN galaxies are known to exhibit bright X-rays, with a mean X-ray luminosity of ∼ 10 43 erg s −1 (e.g., Ricci et al. 2017) similar to that detected by eROSITA.Furthermore, Ricci et al. (2017) and references therewithin showed that a significant fraction of un-obscured AGN, and particularly those selected in optical, tend to exhibit excess soft X-ray emission that can be best described by an absorbed blackbody.They found that this excess can be well fit with a rest-frame blackbody temperature ranging between ∼ 0.5-0.25 keV, with a mean temperature of ∼ 0.1 keV, which is consistent with the blackbody temperature derived by Reusch et al. (2021).Due to the angular resolution of eROSITA, further high resolution X-ray observations would be necessary to confirm whether the detected X-ray emission arises from its host galaxy's AGN or is consistent with the location of AT2019fdr.
If the latter was the case, a detection of X-rays from a SLSN at late times would not be surprising.The total luminosity of the shock and the pre-shock column density of the CSM medium determine the observation features of highenergy radiation.Unless we are in the presence of extremely high shock temperatures or a high ratio of the shock luminosity to the column density, which would guarantee the CSM ionization to a large extent, the photoelectric absorption is an important energy loss channel for high energy photons.The latter could explain the non observation of X-rays at earlier times (Pan et al. 2013).Unfortunately, as already discussed, there could be degeneracies among the parameters, leading to similar properties of the SLSN light curve.Nevertheless, the slow rise of the optical light curve, the very high luminosity peak, and the non observation of X-rays for a considerable amount of time after the first detection would point towards scenarios with highly energetic and relatively low mass ejecta moving through extended high CSM mass stellar winds or shells.

CONCLUSIONS
The IceCube neutrino event IC200530 has been proposed to be in likely coincidence with the source AT2019fdr located at z = 0.2666, observed in the ultraviolet and optical bands, and interpreted as a tidal distruption event candidate in a Narrow-Line Seyfert 1 galaxy.In this paper, we show that the spectra, light curve and color evolution of AT2019fdr may be compatible with the ones of a hydrogen rich superluminous supernova instead.
Under this assumption, the neutrino event IC200530, detected ∼ 300 days after the peak of the electromagnetic emission and with a reconstructed energy of 80 TeV, may have originated as a result of inelastic proton-proton collisions due to the interaction of the supernova ejecta with the circumstellar medium.We find that approximately 10 −8 -5 × 10 −2 muon neutrino and antineutrino events could have been produced by AT2019fdr within the timeframe of interest ( see Table 2), by taking into account the uncertainties on the total ejecta energetics, ejecta mass and on the properties of the the circumstellar medium, as well as the uncertainties in the discrimination of the atmospheric and astrophysical neutrino fluxes.By considering the Eddington bias on neutrino observations, our findings may be compatible with the detection of one neutrino event from AT2019fdr.
In conclusion, the neutrino event IC200530 may be associated with the hydrogen rich superluminous supernova AT2019fdr.As a deeper understanding of the electromagnetic data will become available, neutrinos could be powerful messengers to help to disentangle the nature of AT2019fdr and discover its physics.all the configurations with R CSM = 2 × 10 16 cm and Ẽk = 10 53 erg.As R CSM increases (see the right panel of Fig. 9), the most stringent constraint comes from the compatibility of t diff with the observed light curve.The same trend holds for the case with Ẽk = 5 × 10 52 erg (not shown here), with the difference that there are compatible scenarios with our requirements already for R CSM = 2 × 10 16 cm.For this latter case, for fixed M ej , M CSM and R CSM , the shock velocity v sh is lower, allowing for longer times required to cross the CSM.

B. MAXIMUM PROTON ENERGY
In this appendix, we explore the temporal evolution of E p,max for the set of parameters Ẽk , R CSM , M ej and M CSM considered in this work (see Table 1).We provide an idea of the behaviour of E p,max by displaying in Fig. 10 the ratio between its value at the CSM radius R CSM and the breakout radius R bo , for Ẽk = 10 53 erg with R CSM = 3 × 10 16 cm (left panel) and R CSM = 4 × 10 16 cm (right panel).In both cases, the region where E p,max (R CSM )/E p,max (R bo ) < 1 is the one with relatively low values of M ej /M CSM .Here, either R bo > R dec or R bo R dec , meaning that most of the shock evolution occurs in the decelerating phase (see Eq. 5).When this is the case, the acceleration efficiency drops at a faster rate, leading to decreasing E p,max (see Eq.22).
On the other hand, for large M ej /M CSM , R dec > R CSM is satisfied, implying an increase of E p,max .The intermediate regimes [M ej /M CSM ∼ O(1)] are those in which both free expansion and deceleration occur between R bo and R CSM , being the latter shorter compared to the former, and thus leaving the tendency of E p,max (R CSM )/E p,max (R bo ) to increase unaffected.By keeping Ẽk , M ej and M CSM fixed, a larger R CSM leads to a lower CSM density, with longer t pp ; thus, a larger E p,max (R CSM ) is achievable.This effect is more significant than the slight increase of E p,max (R bo ) for larger R CSM .
Finally, lower values of Ẽk do not compromise the overall trend outlined above.The only effect of decreasing the energy, whilst keeping all other parameters fixed, is to reduce v sh (see Eq. 2) and in turn the acceleration rate, which result in overall smaller values of E p,max .

Figure 2 .
Figure 2. Schematic representation of AT2019fdr after the explosion, assuming spherical symmetry.The central compact object (in black) is surrounded by the SN ejecta (orange region, with the bordeaux arrows indicating the propagation of the ejected material) and a dense CSM envelope (yellow region) which extends up to its outer edge marked by R CSM .The color gradient describes the density gradient (from darker to lighter hues as the density decreases).The dashed black line marks the position of the breakout radius (R bo ).The indigo line represents the forward shock that propagates radially outwards.The black dotted line marks the location of the deceleration radius of the ejecta (R dec ).The latter is located at radii smaller than R CSM (as in this sketch) for a relatively large CSM mass compared to the ejecta mass or larger than R CSM for very massive ejecta and rarefied CSM; see Eq. 4. For extremely large M CSM /M ej , it is possible that R dec < R bo .

Figure 4 .
Figure 4. Muon neutrino and antineutrino fluence from AT2019fdr as a function of the neutrino energy.The reconstructed neutrino energy (E ν ∼ 80 TeV) for IC200530 is marked by a black dotted vertical line.The band encloses the uncertainties on the parameters characterizing AT2019fdr, see Table1.In the proximity of the energy of interest for the interpretation of IC200530, the fluence can vary up to a factor O(10 5 ) in magnitude.Within the allowed parameter space, the lowest fluence is foreseen for configurations with large R CSM , low M CSM and high M ej .The largest neutrino fluence is instead obtained for intermediate values of M CSM and low M ej , which moreover allow a higher proton energy cutoff.

Figure 5 .
Figure5.Muon neutrino and antineutrino event rate expected at the IceCube Neutrino Observatory from AT2019fdr as a function of the time after the shock breakout.The band marks the uncertainty on the neutrino event rate due to the SLSN model parameters, see Table1.The event rate increases rapidly at early times.After peak, the event rates for the SLSN scenarios representing the edges of the envelope decline because of the dominant decreasing trend of v sh as a function of time.In some intermediate scenarios, the increasing trend of E p,max and shallow decrease of v sh can be compensated, providing an increasing event rate at the moment of the detection.The neutrino event IC200530 has been observed ∼ 394 days after t bo as indicated by the dotted vertical line.In the proximity of the detection day, the event rate can vary up to a factor O(10 3 ) in magnitude.

Figure 8 .
Figure8.Contour plot of the total number of muon neutrino and antineutrino events expected at the IceCube Neutrino Observatory from AT2019fdr in 394 days and in the plane spanned by M ej and M CSM for Ẽk = 10 53 erg and R CSM = 4 × 10 16 cm.The black solid lines marks the allowed region of the parameter space, defined by requiring that the location of the shock at the day of neutrino production is still in the CSM envelope and that the SLSN model parameters are compatible with the the light curve rise time.For fixed M CSM , the total neutrino number decreases as M ej increases, given the decline of the shock velocity v sh .Viceversa, for fixed M ej , the number increases as M CSM increases, given the larger number of proton targets for pp interactions.In the region R bo > R dec , one can see an inverted trend.The dotted lines correspond to the contour lines which track the scenarios providing the number of neutrino events displayed, and show how the neutrino number is affected in the transtition from R bo > R dec to R bo < R dec regions.See the main text for more details.

Figure 9 .
Figure9.Left panel: Contour plot of the time the shock takes to travel from R bo to R CSM in the plane spanned by M rj and M CSM .The solid bordeaux line constrains the allowed parameter space by requiring that t CSM − t bo ≥ 394 days (solid bordeaux line).The dashed pink lines constrain the allowed parameter space by requiring that the rising time of the AT2019fdr lightcurve is compatible within a 50% uncertainty with the analytic estimate of the diffusion time provided in Eq. 13; the latter is represented by the solid pink line.Right panel: The same as in the left panel, but for R CSM = 4 × 10 16 cm.For larger R CSM , the crossing time constraint becomes looser, whilst the one related to t diff slowly becomes more stringent.

Figure 10 .
Figure10.Left panel: Contour plot of the ratio between the maximum proton energy E p,max at R CSM = 3 × 10 16 cm and at the breakout radius R bo in the plane spanned by M ej and M CSM .For relatively low values of M ej with respect to M CSM , this ratio tends to decrease.This is due to the fact that for very low M ej /M CSM , R dec < R bo , causing a fast drop of E p,max .Viceversa, for very large M ej /M CSM , the deceleration always occurs at R > R CSM , allowing for a continual increase of E p,max as the time goes by.Intermediate values of M ej /M CSM lead to intermediate trends, with the free expansion and decelerating phase both being present between R bo and R CSM .The dotted black lines indicate the regions the ratio is larger than 1 and 3. Right panel: The same as in the left panel, but for a larger R CSM .The effect of increasing R CSM , while keeping fixed all the other parameters, is to decrease the CSM density and thus to allow for larger E p,max , since the pp interactions become less efficient.

Table 1 .
Benchmark values for the parameters characteristic of AT2019fdr.For the most uncertain ones, we consider a range of variability. 1 bo , solid lines) and at the outer edge R CSM (dashed lines) as functions of the proton energy in rest frame for the SLSN configuration with ( Ẽk , R CSM , M ej , M CSM ) = (10 53 erg, 4 × 10 16 cm, 6 M , 49 M ).The acceleration timescale, pp and adiabatic cooling timescales are represented in red, green and light blue, respectively.Protons are strongly cooled by pp energy losses for all the SLSN parameter configurations considered in this work.