Establishing accretion flares from supermassive black holes as a source of high-energy neutrinos

The origin of cosmic high-energy neutrinos remains largely unexplained. For high-energy neutrino alerts from IceCube, a coincidence with time-variable emission has been seen for three different types of accreting black holes: (1) a gamma-ray flare from a blazar (TXS 0506+056), (2) an optical transient following a stellar tidal disruption event (TDE, AT2019dsg), and (3) an optical outburst from an active galactic nucleus (AGN, AT2019fdr). For the latter two sources, infrared follow-up observations revealed a powerful reverberation signal due to dust heated by the flare. This discovery motivates a systematic study of neutrino emission from all supermassive black hole with similar dust echoes. Because dust reprocessing is agnostic to the origin of the outburst, our work unifies TDEs and high-amplitude flares from AGN into a population that we dub accretion flares. Besides the two known events, we uncover a third flare that is coincident with a PeV-scale neutrino (AT2019aalc). Based solely on the optical and infrared properties, we estimate a significance of 3.6 𝜎 for this association of high-energy neutrinos with three accretion flares. Our results imply that at least ∼ 10% of the IceCube high-energy neutrino alerts could be due to accretion flares. This is surprising because the sum of the fluence of these flares is at least three orders of magnitude lower compared to the total fluence of normal AGN. It thus appears that the efficiency of high-energy neutrino production in accretion flares is increased compared to non-flaring AGN. We speculate that this can be explained by the high Eddington ratio of the flares.


INTRODUCTION
Accreting black holes have long been suggested as potential sources of high-energy particles (Farrar & Gruzinov 2009; Gaisser & Karle ★ E-mail: sjoert@strw.leidenuniv.nl2017) and this expectation was supported by the detection of a highenergy neutrino coincident (at the 3-level) with gamma-ray flaring from the blazar TXS 0506+056 (IceCube Collaboration et al. 2018) and from the nearby active galactic nucleus (AGN) NGC 1068 at the 4-level (IceCube Collaboration et al. 2022).There is also evidence to support neutrino emission from the broader blazar population (see e.g Giommi et al. 2020;Plavin et al. 2020;Hovatta et al. 2021;Kun et al. 2021;Buson et al. 2022), but blazars alone cannot account for the observed high-energy neutrino flux (Aartsen et al. 2017b;Aartsen et al. 2020a;Hooper et al. 2019;Luo & Zhang 2020).Similar to the electromagnetic sky, we expect that the observed cosmic neutrino flux (IceCube Collaboration 2013) arises from multiple source populations (Bartos et al. 2021).
Optical follow-up observations of IceCub neutrino alerts (Aartsen et al. 2017a; Abbasi et al. 2023) using the Zwicky Transient Facility (ZTF; Bellm et al. 2019;Graham et al. 2019;Stein et al. 2023) have identified two optical flares from the centers of galaxies coincident with PeV-scale neutrinos: AT2019dsg (Stein et al. 2021) and AT2019fdr (Reusch et al. 2022).The former belongs to the class of spectroscopically-classified tidal disruption events (TDEs) from quiescent black holes, while the latter originated from a type 1 (i.e., unobscured) AGN (though see Pitik et al. 2022 for an alternative interpretation).The distinctive shared properties we present below suggest these two flares could share a common origin.
Both events show a large amplitude optical flare with a rapid rise time, signalling a sudden increase of mass accretion rate onto the supermassive black hole.Of the ∼ 10 4 AGN detected by ZTF (van Velzen et al. 2021b), less than 1% show similarly rapid and large outbursts (Reusch et al. 2022).The most important unifying signature of the two neutrino-coincident ZTF sources is delayed transient infrared emission, detected by NEOWISE (Wright et al. 2010;Mainzer et al. 2014).This infrared emission is due to reprocessing of the optical to X-ray output of the flare by hot dust ( ∼ 2 × 10 3 K) at distances of 0.1-1 pc from the black hole (Lu et al. 2016;van Velzen et al. 2021a).
A dust reverberation signal is largely agnostic to the origin of the flare near the black hole.Any transient emission at optical, UV or X-ray wavelengths that evolves on a timescale that is shorter than the light travel time to the dust sublimation radius (  ∼ 0.1 pc) will yield a similar-looking dust echo: a flat-topped light curve with a duration of 2  /.This implies that infrared observations of these echoes can be used to construct a sample that unifies "classical TDE" (such as AT2019dsg) and extreme AGN flares (such as AT2019fdr).In this work, we collect a sample of dust echoes from nuclear flares and investigate the significance of their correlation with high-energy neutrinos 1 .
This paper is organized as follows.In section 2 we present the details of the two samples: extreme nuclear transients from supermassive black holes (i.e., accretion flares) and high-energy neutrinos from IceCube.In section 3 we then compute the significance of a correlation between the two samples.This statistical analysis is based only on optical and infrared data.In section 4 we include information from other wavelengths: radio, X-ray, and gamma-rays.In section 5 we discuss the implications of the results.

CATALOG CONSTRUCTION
To build our sample of accretion flares, we select transients from the centers of galaxies as measured with ZTF data and then search for a significant infrared flux increase after the peak of the optical flare using NEOWISE observations.
Below we first present the details of the flare selection and our estimates of the black hole mass.We then present the properties of The lack of ZTF events with a peak after mid-2020 happens because the NEOWISE data release includes observations up to the end of 2020 and to be able to measure the properties of the dust echo we require at least two post-peak detections in NEOWISE.Given the 6 month cadence of NEOWISE, the ZTF light curve has to peak prior to mid-2020 to meet this requirement.Because we allow the neutrino to arrive after the peak of the optical flare, all IceCube events up to the end of 2020 are included in the analysis.
the IceCube neutrino sample and our definition for a flare-neutrino coincidence.

ZTF nuclear flare selection
As described in van Velzen et al. (2019van Velzen et al. ( , 2021b) ) processing of the ZTF alert stream (Masci et al. 2019;Patterson et al. 2019) to yield a sample of nuclear transients is done with AMPEL (Nordin et al. 2019).The input streams include both public ZTF data (MSIP) and private partnership data.We remove events with a weighted hostflare offset (van Velzen et al. 2019) > 0.5".To be able to measure the properties of the light curve we require at least 10 ZTF detections.We also remove ZTF sources for which the majority of the light curve measurements have a negative flux relative to the reference image.These requirements leaves 3142 nuclear transients, see Fig. 1.
To measure the peak flux of the ZTF light curve, we use the observation with the highest flux, after restricting to 90% of the data points with the highest signal-to-noise ratio (excluding 10% of lower quality data makes the peak estimate robust against outliers that occasionally occur in ZTF data).
To define our sample of accretion flares we use a requirement on the rise-timescale and fade-timescale of the ZTF light curve, see Fig. 2.These timescales are obtained from the measurements of van Velzen et al. (2021b) who applied a Gaussian-rise exponentialdecay model to the ZTF alert photometry (both the -and -band are used in this fit).This model explicitly assumes that a single transient explains the entire ZTF light curve.When a light curve has multiple peaks of similar amplitude, the parameters of the fit reflect the (slower) timescale of the majority of the data points.To remove regular variability from normal AGN, we require a minimal amplitude of the flare of Δ < −1 ( denotes the magnitude) and also we set an upper limit to the rise and fade timescale (e-folding time < 150 days and < 500 days, respectively).These cuts leave 1732 sources.The cuts on amplitude and rise/fade times are designed to cast a wide net, recovering all ZTF TDEs and all large-amplitude flares from Seyfert galaxies that have been reported in earlier work The label 'TDE?' indicates accretion flares that occurred in active galaxies (i.e., sources with evidence for accretion prior to the main flare).The three events coincident with a high-energy neutrino are indicated with solid symbols.(van Velzen et al. 2021b;Frederick et al. 2021;Hammerstein et al. 2021).About 15% of the sources are spectroscopically confirmed SN.
To keep track of the follow-up resources and spectroscopic classifications we used the GROWTH Marshal (Kasliwal et al. 2019).Most of the supernova classifications (e.g. as shown in Fig. 2) are based on SEDM (Blagorodnova et al. 2018;Rigault et al. 2019) data obtained for the ZTF Bright Transient Survey project (Fremling et al. 2020;Perley et al. 2020).

NEOWISE dust echo selection
The NEOWISE (Mainzer et al. 2014) light curves cover a period from 2016 to 2020.Most parts of the sky are visited every 6 months and receive about 10 observations within a 24 hour period of this visit (Wright et al. 2010).The inverse-variance weighted mean of the cataloged flux during these visits is used to construct the NEO-WISE light curves.For each source, the baseline is defined using all NEOWISE observations obtained up to 6 months before the peak of the ZTF light curve.The 6 months padding is added to avoid including part of the dust echo signal into the baseline (e.g., when ZTF observations miss the onset of the flare).We see that the three accretion flares coincident with a high-energy neutrino (filled symbols) are among the strongest dust echoes in the sample of nuclear transients.
To measure the echo strength we require two observations after the ZTF peak.We define the dust echo flux as Δ IR : the difference between the baseline flux and the mean NEOWISE flux within one year after the ZTF peak.The echo strength is Δ IR / rms , where  rms is the root-mean-square variability of the baseline observations.The significance of the rms variability is measured using the ratio  rms /  , where   is the measurement uncertainty of the baseline observations.We selected candidate dust echoes by requiring that the echo strength is larger than the significance of the baseline rms variability: Δ IR / rms >  rms /  .We apply this criterion to the light curves of both NEOWISE bands (W1 and W2; central wavelengths of 3.4 and 4.6 m, respectively).This selection leaves 140 nuclear transients with candidate dust echoes.After selecting accretion flares based on the ZTF properties (Fig. 2), we are left with 63 flares with candidate dust echoes (Table 3).In Fig. 3 we show the echo flux and luminosity (for the subset of sources with a spectroscopic redshift) versus echo strength.
The time difference between the optical and infrared light curve peaks yields an estimate of the inner radius of the dust reprocessing region.At this dust sublimation radius, the bolometric flux absorbed by the dust is equal to the infrared luminosity emitted by the dust (with a spectrum that is determined by the sublimation temperature of the dust).We can therefore estimate the bolometric luminosity of the flare from the duration of the infrared reverberation light curve (Lu et al. 2016;van Velzen et al. 2016van Velzen et al. , 2021a)).Using this geometric luminosity estimate (Eq.12 in van Velzen et al. 2016), we find a bolometric luminosity ∼ 10 45 erg s −1 for all accretion flares coincident with high-energy neutrinos.All are consistent with reaching the Eddington limit (Table 1).For infrared emission due to reverberation, the energy emitted by the dust cannot exceed the integrated bolometric luminosity of the flare.For this reason, the lower optical-to-infrared ratio of the third source (cf.Fig. 13) likely implies a larger bolometric correction for its optical emission.This suggests ≈ 1 mag of optical extinction.

Black hole mass estimates
For all nuclear flares in our sample, black hole masses are can be estimate if optical spectra are available.We either use a relation based on reverberation mapping (Blandford & McKee 1982;Peterson et al. 2004), also known as "virial mass estimates", or the - relation (Magorrian et al. 1998;Gebhardt et al. 2000).We use the former for spectra of type 1 AGN, i.e., sources that show broad Balmer emission lines in their optical spectrum.The - relation is used for sources without broad emission lines that have a host galaxy spectrum with a measurement of the velocity dispersion of the stars.
For the reverberation method, a measurement of the size of the broad-line region in combination with the velocity of the emission lines yields an estimate of the black hole mass.This distance of the broad-line region to the black hole is not measured directly, but follows from the observed disk luminosity.We adopt the relation from Ho & Kim (2015): 5100 10 44 erg s −1 0.533 + 6.91 (1) Here  5100 is the continuum luminosity at 5100 Å in the rest frame.For active galaxies with spectra from the Sloan Digital Sky Survey (SDSS; York et al. 2000) we use the catalog of Liu et al. (2019), who selected 14,584 type 1 AGN based on detection of a broad H line and applied Eq. 1 to measure the black hole mass.We find 580 nuclear transients with black hole mass estimates based on this catalog.Of these, eight are classified as accretion flares with potential dust echoes.In addition, spectroscopic follow-up observations of ZTF transients have yielded nine more type 1 AGN, six of which have been published (Frederick et al. 2019(Frederick et al. , 2021) ) and three are presented for the first time in this work: AT2019meh (ZTF19abclykm), AT2020afab (ZTF19abkdlkl), and AT2020iq (ZTF20aabcemq).We also obtained a new post-peak spectrum of AT2019aalc (Fig. 4), which shows ongoing accretion two years after the peak of the optical flare.
For the remaining transients that have a velocity dispersion measurement based on a spectrum of the host galaxy, we apply the Gültekin et al. (2009) - relation.This yields 219 additional  BH measurements, of which five are classified as accretion flares with potential dust echoes.Of these five, three are based on archival SDSS spectra of the host galaxy (Thomas et al. 2013) and two are based on follow-up observations obtained after the flare has faded (Nicholl et al. 2020;Cannizzaro et al. 2021).In Table 3, we list the reference for the black hole mass estimate for each accretion flare.

IceCube neutrino alerts
Our parent sample of neutrino alerts (Aartsen et al. 2017a) includes the events published up to December 31, 2020.We exclude the Ice-Cube alerts that were subsequently retracted after the automated alert was issued.We also remove two events without a reported signalness (IC190331A; Kopper 2019 and IC200107A; Stein 2020a) and one event with a 90%CL area larger than 300 deg 2 (IC200410A; Stein 2020c).This leaves 43 events, listed in Table 2. Of these, three fall outside ZTF extra-galactic footprint (Ω ZTF = 2.8 × 10 4 deg 2 ; Stein et al. 2021), leaving 40 IceCube alerts that can yield a coincidence with an accretion flare in our ZTF+NEOWISE dataset.
The signalness, or  astro , measures the probability that a detector track recovered by IceCube is from a cosmic neutrino, based on the reconstructed energy and assumptions about the astrophysical neutrino flux.The sum of signalness of the 40 neutrinos in the ZTF footprint is 17.7.Multiplying by 0.9 to account for the 10% of neutrinos whose true location is outside the reported 90%CL area, we estimate that about 16 cosmic neutrinos can in principle be recovered by our analysis.

Flare-neutrino coincidence
An accretion flare is considered coincident with a neutrino when the source falls inside the 90%CL reconstructed neutrino sky location and this flare is detected in ZTF and NEOWISE when the neutrino arrives, with a maximum delay of one year relative the peak of the optical light curve.For longer delays, our search would lose sensitivity because our NEOWISE dataset only contains photometry up to the end of 2020 (see Fig. 1).This requirement for coincidence yields three matches between the 43 IceCube alerts (section 2.4) and the 63 accretion flares with potential dust echoes (section 2.2).
We immediately notice that the three accretion flares with a coincident neutrino have very strong dust dust echoes compared to the rest of the nuclear flare population (5).The significance of this result is discussed below in section 3) Black hole mass estimates based on optical spectroscopy are available for about 1/3 of the accretion flares (see section 2.3).We find that strong dust echoes (Δ IR / rms > 10) are observed almost exclusively for events with  BH < 10 8  ⊙ , see Fig. 6.The threshold for strong echoes appears to coincide with the maximum mass for a visible disruption from a solar-type star (Hills 1975).This motivates the construction of a TDE candidate sample: all accretion flares with  BH < 10 8  ⊙ .The infrared properties of the accretion flares with black hole mass estimates are used to place all nuclear transients with dust echoes into a two-tier classification scheme: (1) accretion flares with strong dust echoes and (2) regular AGN variability.In the next section, we consider the hypothesis that the former class is a source of high-energy neutrinos.

STATISTICAL SIGNIFICANCE
We measure the strength of the observed flare-neutrino correlation by defining a test statistic based on the likelihood ratio.Our statistical test accounts for both the spatial localization of the neutrino and the infrared properties of the flare relative to the TDE candidate  population.To quantify the statistical significance of our result we randomly redistribute the accretion flares into the ZTF footprint and compute the test statistic for a large number of these simulations.Similar to the approach used for the neutrino detection from the blazar TXS 0506+056 (IceCube Collaboration et al. 2018), our likelihood analysis is not blind because the input was defined after two neutrino associations had been established (AT2019dsg and AT2019fdr).However, we note that these prior associations were based on ZTF follow-up observations of neutrino alerts.The infrared data relevant for our analysis were released in March 2021, and therefore the echo strength was not used to establish the neutrino coincidence.Below we present the details of our likelihood ratio test (section 3.1), including its statistical power (section 3.2) and robustness (section 3.3) we also demonstrate how the significance can be estimated without the use of Monte Carlo simulations (section 3.4).The software and the data required to reproduce these likelihood estimates can be obtained via Zenodo (van Velzen & Stein 2022).

Likelihood-ratio based estimate of significance
To test for a correlation between high-energy neutrinos and accretion flares we consider the likelihood ratio between a signal hypothesis and a background hypothesis.For a given neutrino, , there are are several possible hypotheses that can be used to explain the origin.If there are  accretion flares, we can test  discrete hypotheses, i.e., that the neutrino originated from source #1 ( 1 ), from source #2 ( 2 ), etc.We can denote this group of hypotheses as signal hypotheses, i.e.,   ≡   .Additionally, we have a further possible explanation, namely that the neutrino did not originate from any of the accretion flares.The neutrino itself might not be of astrophysical origin (but instead is due to the atmospheric background), or the neutrino may originate from an astrophysical source that is not included in our sample of accretion flares.All these alternative options are included in the hypothesis  0 .We can also denote this null hypothesis, or background hypothesis, as  0 ≡ .
For each hypothesis, we can define the likelihood that the data, , is obtained under that hypothesis L  = (|  ).For each neutrino, we select the hypothesis with the greatest likelihood as our best fit: In the rare occasion that multiple accretion flares are found in coincidence with a single neutrino, we select the flare with the highest likelihood: Next, we compare the signal hypothesis to the background hypothesis: When no accretion flare is found in coincidence with a neutrino, we have L = L 0 .Below, we first discuss the input for the signal hypothesis, followed by the components of the background hypothesis.

Components of the signal hypothesis
The probability of the data under a particular signal hypothesis is given by the joint probability of two distinct components: (i) The probability that a spatial and temporal coincidence between a given neutrino alert and accretion flare would be observed if the neutrino was produced by the accretion flare.We denote this with  coin (, , |), with  and  the right ascension and declination of the accretion flare, respectively.
(ii) The probability that the properties of the accretion flare would be observed, if the neutrino was produced by the accretion flare.
Probability 1 of this list is trivial to calculate, because spatial coincidence is established when the accretion flare falls inside the reported 90%CL localization area of the neutrino.Hence for coincident events,  coin = 0.9.Here we implicitly assume that all neutrinos emitted by accretion flares arrive within our temporal search window.
Probability 2 of the list yields a second and third term that enter the likelihood of the signal hypothesis.These are based on the echo strength and the echo flux.We motivate these two terms below.
A key property of our signal hypothesis is that stronger echoes are less likely to be explained by normal AGN variability (Fig. 6), but rather by a single outburst or transient.We estimate the PDF of the echo strength from the observed distribution of the 18 accretion flares that are candidate TDEs (i.e., sources with an estimate black hole mass < 10 8  ⊙ ).To turn the binned distribution into a continuous PDF, we use a Gaussian kernel density estimate (KDE).The result is shown in Fig. 7.For all KDE estimates, we select the optimal bandwidth following Scott's Rule (Scott 2015).
The For the signal hypothesis we assume  (ΔF IR |S) ∝ ΔF IR (purple line).The observed echo flux of the three events that found in coincidence with an IceCube neutrino alert is consistent with this PDF (See Fig. 12).
production: a linear coupling between the total electromagnetic luminosity and the neutrino luminosity.To estimate the neutrino flux at Earth, we thus need the total electromagnetic fluence (i.e., the bolometric energy over the distance squared).The bolometric energy cannot be measured directly because a large (but unknown) fraction of the energy is emitted at higher frequencies than the optical wavelength range of ZTF.Fortunately, dust absorption is efficient up to soft X-ray frequencies, hence the dust reprocessing light curve provides a good tracer of the bolometric luminosity.With these simplifying assumptions, we thus obtain a linear scaling of the neutrino flux with infrared flux of the echo: Here the max/min correspond to the maximum/minimum observed echo flux of the accretion flares, respectively.After applying our test statistic, we will validate the consistency of this linear scaling.To conclude, the signal likelihood of a single neutrino and accretion flare coincidence is given by:

Components of the background hypothesis
We now consider the probability of observing the data under the null hypothesis.We again have two distinct factors: (i) The probability of a spatial and temporal coincidence if the neutrino was not produced by any accretion flare:  coin (, , |).
(ii) The probability that the properties of an accretion flare found in coincidence would be observed if the neutrino was not produced by that accretion flare.
The first item in this list accounts for the so-called chance coinci-dences, both spatially and temporally.The expectation value for the number of chance coincidences within the 90%CL area (Ω  ) of a given neutrino is given by  bg × Ω  , with  bg being the areal density of coincident events that are expected for the background hypothesis.Because this expectation value for a single neutrino is always ≪ 1, the Poisson probability of obtaining a single spatial and temporal chance coincidence between a neutrino and an accretion flare can be written as A Monte Carlo approach is used to estimate  bg .We need to assign a peak time and location for each simulated accretion flare.Since our window for temporal coincidence is broad, we can use a non-parametric method to assign a time of peak, namely shuffling the peak times for each Monte Carlo realization.Next, we need a method to simulate celestial coordinates of ZTF extra-galactic tran-sients.We start from the parent sample of 3142 nuclear transients with at least two post-peak NEOWISE observations (see section 2.1).The areal density in this sample is too low to allow a non-parametric approach.We therefore applied Gaussian KDE to obtain a continuous two-dimensional PDF of the celestial coordinates.To ensure the lack of events from the Galactic Plane is properly reflected, we enforce zero probability for Galactic coordinates with || < 8 deg.After this small modification, we can use the resulting PDF to simulate celestial coordinates of ZTF nuclear transients using rejection sampling.These steps are summarized in Fig. 9.After simulating 10 6 samples of 63 accretion flares, we find a total of 4.7 × 10 5 coincident neutrinos.Hence the total number of expected coincident events in a sample of 63 accretion flares is 0.48.To obtain the areal density of background coincidences, we divide this expectation value by ΣΩ  = 698.6 deg 2 , the total neutrino area that overlaps the ZTF footprint:  bg = 0.48/ΣΩ  = 6.9 × 10 −4 deg −2 .
The second probability in the list above can be found by considering the properties of accretion flares expected for a chance coincidence.Such accretion flares will be drawn at random from the general population.The probability to detect a given dust echo flux can be estimated from the flux distribution of the parent population of all nuclear transients with detected transient infrared emission.The probability to detect a given echo strength follows from the distribution of nuclear transients with a potential echo.Because the TDE candidate population dominates the upper-end of the echo strength distribution, we exclude these flares when applying a Gaussian KDE to obtain (Δ IR / rms |).The resulting PDFs are displayed in Figs.7 & 8.

Final result: weighted likelihood ratio test
In principle, the likelihood ratio test could be extended to include the probability of observing the neutrino properties detected at the IceCube Neutrino observatory for the background and signal hypothesis.This would account for the fact that some detector properties, such as detected energy, can be used to infer the probability of a neutrino being astrophysical.However, this information is not provided by IceCube.Instead, IceCube provides an estimated astrophysical probability for each neutrino event based on assumed properties of the astrophysical neutrino flux ( astro , see section 2.4).We therefore choose to simply use this astrophysical probability as a weight in our likelihood ratio test.
We now collect the signal and background terms to define our test statistic: The sum runs over all neutrinos and the flare properties are evaluated for the (best-fit) accretion flare that is found in spatial and temporal coincidence with the neutrino.
For the observed sample of 63 accretion flares we find TS= 29.6.We can now estimate the significance of this result by computing the TS distribution under the background hypothesis.We use 10 6 samples of 63 accretion flares drawn from the PDF of the celestial coordinates of nuclear transients (see Fig. 9), with shuffled flare start times.For each of these Monte Carlo realizations, we compute the TS using the coincident neutrino and flare pairs.The resulting distribution of TS is shown in Fig. 10.A fraction 1.5 × 10 −4 of the simulations for the background hypothesis have an equal or greater TS than the observations, corresponding to a significance of 3.6.The detection of a third neutrino in coincidence with an accretion flare (i.e., AT2019aalc), decreases the probability of the background hypothesis by a factor 60; if our search had not uncovered this new event, the significance would have been 2.4.

Statistical power
It can be instructive to investigate the statistical power of our likelihood method (Eq.8) as a function of the expectation value of the number of signal neutrinos,  exp .To simulate a signal for a given  exp , we draw  events from a Poisson distribution with expectation value  exp .Each of these  events has a 90% probability to be detected in coincidence with an accretion flare.To each simulated neutrino-flare pair we assign an echo flux and strength from the PDF of the signal hypothesis (see Figs. 7 & 8).We repeat this signal simulation 10,000 times to obtain the TS distribution for a given  exp .We compare this distribution to a few critical TS values: TS = 0, which is the median of the background TS distribution, and TS = 27.9,TS = 43.8, and TS = 57.5 which correspond to the threshold for a 3, 4, and 5 significance, respectively (the limiting TS value for 5 is obtained by extrapolating the background TS distributing, using the approach outlined in Aartsen et al. 2017c).
When at least 50% of the simulated signal trials have TS > 0, our test statistic can be called admissible (on average, we have enough sensitivity to prefer the signal hypothesis over the null hypothesis); this happens when  exp > 0.3.For reference, our detection of three coincident neutrinos implies  exp = 3 +3.8−1.9 (90% CL).From Fig. 11 we see that for this range of  exp , a significance between 3 and 4 should be expected.
The number of signal neutrinos is related to the fraction of cosmic high-energy neutrinos that originate from accretion flares,  flare : Here  cosmic ≈ 16 is the number of astrophysical neutrino alerts Fraction of cosmic neutrinos that orignate from accretion flares Figure 11.Assessment of the statistical power of our weighted likelihood ratio test (Eq.8).The horizontal axis encodes the expectation value for the number of neutrinos that are detected in coincidence with an accretion flare in our sample ( exp in Eq. 9).For each expectation value, one can simulate a distribution of the test statistic (TS; Eq. 8) and compute the fraction of injected signal trials above a given TS threshold value.The top axis shows the fraction of astrophysical IceCube alerts that originate from accretion flares in our catalog (based on Eq.9).For reference, based on the detection of three coincident neutrino-flare pairs, the 90%CL range of  exp is (1.1, 6.8).
that are included in our search (see Section 2.4).The parameter  ZTF accounts for the fraction of astrophysical neutrinos from accretion flares that detected by IceCube, but not included in our catalog.In particular, a large population of relatively high-redshift accretion flares ( > 0.5) will not be detected in ZTF or NEOWISE, but these can yield a sizeable fraction of neutrino alerts.Following Stein et al. (2021), we adopt  ZTF = 0.5.The resulting  ac is shown at the top of Fig. 11.

Cross-checks and robustness
Now that we have established evidence for three accretion flares with a neutrino counterpart, we can do a cross-check on the assumption that the neutrino flux is coupled linearly to the infrared echo flux.We first make an ansatz for the coupling strength between the neutrino expectation value (  ) and the infrared flux:   = Δ IR , with  = 100Jy −1 .Applying this coupling to the observed echo flux distribution of the TDE candidates we obtain   for each TDE candidate.
After simulating neutrino detections from each TDE by drawing from a Poisson distribution, we find that, for this coupling strength, 36% of the simulations yield at least three difference sources, each with at least one neutrino.In most cases, we obtain one neutrino per source (e.g., for the subset of simulations with exactly three neutrino sources, 46% yield more than one neutrino from a single flare).In Fig. 12 we show the distribution of echo flux for the simulations that yield exactly three detected neutrinos from three different sources.This prediction is consistent with the flux distribution of the three flares that have evidence for a neutrino counterpart.While this agreement is encouraging, we should expect deviations from some of the simplifying assumptions that are made to construct the test statistic (Eq.8).This will decrease the sensitivity of our statistical method.Below we test the robustness of our significance estimate by making several alternative choices for the input of Eq. 8.
First, the simulation of the expected echo flux distribution for a linear coupling with the neutrino expectation value and infrared flux could also be considered as input for (Δ IR |).This approach is not ideal for two reasons: (1) the coupling strength is a free parameter; (2) the numbered of observed coincident events is needed to normalize this PDF.Nevertheless, if we redo our estimate of the significance using the predicted dust echo flux distribution as shown in Fig. 12, we obtain  = 1.7 × 10 −4 for the background hypothesis (3.6).
Next, we consider the PDF of the celestial coordinates of nuclear transients, which is required to generate the distribution of background flares in our Monte Carlo analysis.The resulting significance is not sensitive to the details of this PDF.To demonstrate this, we can draw the background Monte Carlo samples from the observed local galaxy distribution, without KDE smoothing.We use the 2MASS redshift survey (Huchra et al. 2012) and select all galaxies that fall within the declination limits of ZTF (−28 <  < 86), yielding 3.1 × 10 4 sources.The resulting sky distribution is uniform within the ZTF footprint, which is not consistent with the observed distribution of nuclear transients (Fig. 9).While the 2MASS galaxy distribution is clearly not an appropriate description of the celestial coordinates of our nuclear flare sample, this has a modest influence on the inferred significance if this map is selected to generate background samples.After drawing the coordinates of the background samples from the 2MASS galaxies, we find  bg = 8.8 × 10 −4 and Finally, we consider the impact of the selection of accretion flares based on the ZTF properties of the light curve.If we make no cuts on the ZTF properties and simply select all potential dust echoes (i.e., Δ IR / rms >  rms /  ) we obtain 163 sources.For this sample, we expect 1.9 matches under the background hypothesis (compared to 0.48 for the original sample of 63 accretion flares).This larger population yields one additional coincident neutrino (for this fourth source, ZTF19aaozcxx, the background hypothesis is preferred over the signal hypothesis, TS  = 0).The probability of the null hypothesis is  = 3.2 × 10 −4 (3.4).

Density-based estimate of significance
As a final cross-check on the robustness of our likelihood and Monte Carlo method we consider a simple, yet instructive estimate of the significance based solely on the areal density of accretion flares with large dust echoes.For this estimate, no Monte Carlo sampling is needed.We first note that for Δ IR / rms > 10, the flare population is dominated by TDE-like echoes (Fig. 6).Applying this cut on echo strength leaves 29 flares.The fraction of neutrino alerts that yield a temporal coincidence with these flares is  temp = 0.24.We thus obtain the effective source density of large echoes:  eff = 29 ×  temp /Ω ZTF = 2.5 × 10 −4 deg −2 , with Ω ZTF = 2.8 × 10 4 deg 2 the extragalactic sky seen by ZTF (Stein et al. 2021).Multiplying this effective source density with the total area of the IceCube neutrino alerts, we obtain the expected number of chance coincidences.This expectation value is 0.17 and the Poisson probability to see at least three events is 8 × 10 −4 (3.2).

Detection at radio wavelengths
The first neutrino-detected source AT2019dsg, showed transient radio emission, evolving on a timescale of month, with a peak luminosity (  ) of 8 × 10 38 erg s −1 at 10 GHz (Stein et al. 2021).The source AT2019fdr was detected in radio follow-up observations with similar luminosity (2 × 10 39 erg s −1 at 10 GHz), with marginal evidence for variability at 10 GHz frequencies (Reusch et al. 2022).Due to its higher redshift, the radio emission from AT2019fdr is relatively faint (0.1 mJy) and below the detection threshold of wide-field radio surveys such as FIRST (Becker et al. 1995) or VLASS (Lacy et al. 2020).
Finally, the third accretion flare coincident with a high-energy neutrino, AT2019aalc, is detected in both FIRST and VLASS.The FIRST observation were obtained in the year 2000, thus predating the optical flare by almost two decades.These radio observations yield a luminosity of 2 × 10 38 erg s −1 at 1.4 GHz.The VLASS observations were obtained on two dates, 2019-03-14 and 2021-11-06, yielding 3 GHz radio luminosities of 3 × 10 38 and 5 × 10 38 erg s −1 , respectively.This factor ≈ 2 flux increase from the first VLASS epoch (three months before the optical peak) to the second VLASS epoch (2 years post-peak) is statistically significant (at the 8-level, as estimated from the rms in the VLASS "Quick Look" images).

Fermi gamma-ray upper limit
We analysed data from the Fermi Large Area Telescope (Fermi-LAT; Atwood et al. 2009), a pair-conversion telescope sensitive to gamma rays with energies from 20 MeV to greater than 300 GeV.Following the approach outlined in Stein et al. (2021), we use the photon event class from Pass 8 Fermi-LAT data (P8R3_SOURCE), and select a 15 × 15 deg 2 region centered at the target of interest, with photon energies from 100 MeV to 800 GeV.We use the corresponding LAT instrument response functions P8R3_SOURCE_V2 with the recommended spectral models gll_iem_v07.fitsand iso_P8R3_SOURCE_V2_v1.txtfor the Galactic diffuse and isotropic component respectively, as hosted by the FSSC.We perform a likelihood analysis, binned spatially with 0.1 deg resolution and 10 logarithmically-spaced bins per energy decade, using the Fermi-LAT ScienceTools package (fermitools v1.2.23) along with the fermipy package v0.19.0 (Wood et al. 2017).
We studied the region of AT2019aalc in the time interval that includes 207 days of observations from the discovery of the optical emission on 2019 April 26 to the observation of the high-energy neutrino IC191119A on 2019 November 19.In this time interval, there is no significant (≥ 5) detection for any new gamma-ray source identified with a localization consistent with IC191119A.Two sources from the fourth Fermi-LAT point source catalog (4FGL-DR2; Abdollahi et al. 2020;Ballet et al. 2020), consistent with the IC191119A localization, are detected in this interval.These are 4FGL J1512.2+0202(4.1 deg from IC191119A), associated with the object PKS 1509+022, and 4FGL J1505.0+0326(5.0 deg from IC191119A), associated with the object PKS 1502+036.The flux values measured for these two detections are consistent with the average values observed in 4FGL-DR2.
Likewise, we studied Fermi-LAT data of AT2019fdr using a time window from its discovery on 2019 May 3 to the arrival time of IC200530A on 2020 May 30.In this window, no gamma-ray sources were significantly (≥ 5) detected within the localization region of IC200530A, including both previously known 4FGL-DR2 catalogued sources and new gamma-ray excesses.
For both AT2019fdr and AT2019aalc, we test a point-source hypothesis at their position under the assumption of a power-law spectrum.The 95% CL upper limit for the energy flux (i.e., integrated over the whole analysis energy range) listed in Table 1 is derived for a power-law spectrum (/ ∝  −Γ ) with photon power-law index Γ = 2.The Fermi-LAT upper limit for AT2019dsg listed in Table 1 is obtained using the same LAT data analysis strategy and covers a similar time window (150 days) relative to the optical discovery and the neutrino arrival time (Stein et al. 2021).
No significant gamma-ray emission is detected in Fermi-LAT data at the source positions prior to their optical discoveries (Abdollahi et al. 2017).

SRG/eROSITA X-ray detections
The SRG X-ray observatory (Sunyaev et al. 2021) was launched to the halo orbit around Sun-Earth L2 point on 2019 July 13.On 2019 December 8, it started the all-sky survey, which will eventually comprise eight independent 6-month long scans of the entire sky.In the course of the sky survey, each point on the sky is visited with a typical cadence of 6 months (the exposure and number of scans depends on ecliptic latitude).The eROSITA soft X-ray telescope (Predehl et al. 2021) operates in the 0.2-9 keV energy band with its effective area peaking at ≈ 1.5 keV.
As of October 2021, AT2019aalc (= SRGe J152416.7+045118) has been visited four times starting from February 2, 2020 with 6 month intervals and was detected in each scan.The X-ray light curve of the source as seen by eROSITA reached a plateau between 2020 August and 2021 January with the 0.3 − 2.0 keV flux of ≈ 4.6 × 10 −13 erg s −1 cm −2 .The source had a soft thermal spectrum with the best fit blackbody temperature of  = 172 ± 10 eV.
The source AT2019fdr (= SRGe J170906.6+265124) has been visited four times starting 2020 from March 13.The sources was detected only once, on 2021 March 10-11, with a 0.3-2.0keV flux of ≈ 6.0 × 10 −14 erg s −1 cm −2 and an extremely soft thermal spectrum with a temperature of 56 +32 −26 eV.This flare displayed the softest X-ray spectrum of all TDEs detected by eROSITA so far (Sazonov et al. 2021).In the three visits when the source remained undetected, the upper limit on its flux was in the range of ∼ (2 − 5) × 10 −14 erg s −1 cm −2 , see Reusch et al. (2022) for details.
The source AT2019dsg has been visited by eROSITA three times starting from 2020 May 9 and so far remained undetected with the  2).The upper limit on the neutrino luminosity is estimated by assuming an expectation value of one particle in Δ (the true neutrino luminosity could be one or two orders of magnitude lower due to the Eddington bias, see Sec. 5.2).The X-ray temperature (  ) and luminosity (  ) are based on SRG/eROSITA or Swift/XRT and cover an energy range of 0.2 to 10 keV (see Sec. 4.3).The radio luminosity is measured at  ≈ 5 GHz (see Sec. 4.1).The 95% CL upper limit on the -ray luminosity (  ) is obtained from data of the Fermi-LAT telescope and covers photons in the energy range 0.1-800 GeV (see Sec. 4.2).The black hole mass ( BH ) is estimated from the optical spectrum of the host galaxy (Sec.2.3).The Eddington ratio (  Edd ) follows from the bolometric luminosity ( bol ) as estimated from the duration of the dust echo (Sec.2.2).
3 upper limit for the combined data of the three visits of ≲ 1.9 × 10 −14 erg s −1 cm −2 (assuming a power law spectrum with the slope of Γ = 1.8).The X-ray measurements of AT2019dsg listed in Table 1 are based on Swift/XRT; Gehrels et al. 2004) observations that were obtained closer to the optical peak of the flare (van Velzen et al. 2021b;Stein et al. 2021) than the eROSITA observations.

A new population of neutrino sources
By using infrared observations of dust echoes as a tracer of large accretion events near black holes, we are able to unify TDEs in quiescent galaxies and (extreme) AGN flares.This allows us to test a new hypothesis: large amplitude accretion flares are sources of highenergy neutrino emission.Thanks to our systematic selection of dust echoes, we can use a large sample of 40 neutrinos, compared to the 24 neutrinos that were followed-up by ZTF to find AT2019dsg and AT2019fdr (Stein et al. 2021;Reusch et al. 2022).This larger sample allows us to uncover a third flare (ZTF19aaejtoy/AT2019aalc), which happens to have the highest dust echo flux of all ZTF transients.The significance of this population of three flares is 3.6.The detection of the third flare reduces the probability of a chance association by a factor of ≈ 60.Additional evidence for neutrino emission from accretion flares follows from the shared multi-wavelength properties of the three sources with a neutrino counterpart.
For AT2019aalc, we measure a soft thermal spectrum with temperature of  = 172 ± 10 eV.Such soft thermal emission is rare: of all accretion flares with potential dust echoes, less than 13% are as soft as AT2019aalc ( < 172 eV), and even fewer are as soft as AT2019dsg and AT2019fdr.These distinctive X-ray properties provide 3-level ( = 0.13 3 ) evidence for the hypothesis that accretion flares are correlated with high-energy neutrinos.
Another shared property of the three events is low-luminosity radio emission (see section 4.1).In our sample of accretion flares with dust echoes, less than 10% are detected in archival radio observations (FIRST or VLASS) and a similarly low fraction of TDEs is detected in radio follow-up observations (Alexander et al. 2020).If the three neutrino associations in our sample of accretion flares are due to chance, the probability to find three radio detections is < 10 −3 .
Transient radio emission and soft X-ray emission is not consistent with a supernova explanation for AT2019fdr that has been proposed by Pitik et al. (2022).Taken together, the shared X-ray, radio, and optical properties of the three flares with neutrino coincidences (Fig. 13) point to a new population of cosmic particle accelerators powered by transient accretion onto massive black holes.

Energetics
Around PeV energies, the expectation value for the number of Ice-Cube neutrino detections can be approximated as with  neutrino the total energy carried by (mono-energetic) neutrinos (Stein et al. 2021).Our single-neutrino associations suffer from a significant Eddington bias (Strotjohann et al. 2016).Since our analysis considered a sample of 63 accretion flares, we should expect  neutrino < 1 for any individual source in our sample.This implies Eq. 10 cannot be inverted to estimate the energy in neutrinos emitted by a single source.Only a upper limit can be estimated by setting  neutrino = 1, the result is shown in Table 1.
Instead of investigating the neutrino luminosity of a single source, a more useful approach is to demand that the expectation value for the number of neutrino associations for the full population of 63 accretion flares is equal to the observed value (i.e., three).To obtain  neutrino for the entire population, we again assume a single coupling strength (  ) between the electromagnetic energy and the energy carried by neutrinos: We estimate the bolometric energy of the flare from its bolometric luminosity as obtained from the duration of the dust echo ( bol , see Sec. 2.2) and the -folding time  of the decaying part of the optical light curve (see Fig. 2):  bol =  bol .The total bolometric flare energy of the 63 accretion flares is 1 × 10 54 erg.Hence an expectation of  neutrino = 3 +4 −2 (90% CL) yields a mean coupling strength   = 4 +5 −3 × 10 −2 .For most models of particle acceleration in AGN accretion disks, the luminosity in high-energy cosmic rays is between one and two orders of magnitude larger than the neutrino luminosity (Begelman et al. 1990;Murase et al. 2020).Our estimate of the fraction of the electromagnetic energy that is given to neutrinos (  ∼ 10 −2 ) thus implies that the energy in cosmic rays could be of the same order as the bolometric luminosity emitted by the flare.
If we apply our simple model (i.e., a fixed coupling strength between the bolometric energy and the energy emitted in neutrinos) to individual sources, we find a neutrino expectation value of 0.08, 0.004, and 0.21 for AT2019dsg, AT2019fdr, and AT2019aalc, respectively.We note that for AT2019dsg, a model of particle acceleration in the core of an accretion disk with a constant accretion rate yields an expectation of 0.1 IceCube high-energy neutrinos for an integration time of 1 year (Murase et al. 2020).We refer to Winter & Lunardini 2023 for further exploration of potential particle acceleration mechanisms in all three accretion flares.

Neutrino arrival delay
If the optical luminosity and the luminosity of high-energy neutrinos are directly coupled, the fraction of the optical energy that is emitted at the time the neutrino arrives should follow a uniform distribution.Within our search window of 1 yr after the optical peak, this fraction is 0.85, 0.87, and 0.66 for AT2019dsg, AT2019fdr, and AT2019aalc, respectively.A Kolmogorov-Smirnov test yields  = 0.08 for the null hypothesis that these values are drawn from a uniform distribution.
While the relatively late arrival time of the three neutrinos is not statistically significant, the effect is large enough to permit some speculation about possible origins for delayed neutrino emission in TDEs or flaring AGN.Below we suggest three possible explanations.
First, if the mass accretion rate is constant for about one year, the neutrino flux can also be expected to be constant over this period.A constant neutrino flux implies that a delayed neutrino detection is equally likely as a detection close to the optical peak.This is possible because the optical/UV emission of TDEs might not trace the accretion rate (Piran et al. 2015) and instead is emitted at the first shock of the stellar debris streams, which marks the onset of the debris circularization (Bonnerot & Stone 2021), see Roth et al. (2020) for a review.A roughly constant accretion rate following the first year after disruption was used to explain the delayed neutrino detection of AT2019dsg, and this idea is supported by the radio (Stein et al. 2021) and X-ray observations of this event (Mummery 2021).
Second, the mass accretion rate may not be constant, but delayed.The circularization timescale of the stellar debris can be estimated as (Hayasaki & Jonker 2021).Here  represents how efficiently the kinetic energy at the stream-stream collision is dissipated; the most efficient ( = 1) case corresponds to the result of Bonnerot et al. (2017).Because the inner accretion disk is formed after the circularization of the debris, it could take of order  circ for the first neutrinos to be produced, which appears to match the observed time delay (Fig. 13).
A caveat to Eq. 12 is a potential interaction of stellar debris stream with a pre-exisiting accretion disk, which can shorten the circularization timescale (Chan et al. 2019).Another caveat is that the radiative efficiency of the first stream-stream shock might be too low to explain the prompt optical/UV emission (Lu & Bonnerot 2020).Instead, the early-type optical/UV emission could be attributed to reprocessing of photons emitted from the accretion disk, implying the disk has already formed when the optical emission peaks (Bonnerot & Lu 2020).In this case, for PeV particle acceleration in the newly formed accretion disk, the debris circularization time is too short to explain the neutrino arrival delay.
Finally, a delay between optical emission and the neutrino arrival time can be explained if particle acceleration happens in a jet or outflow and neutrinos are only produced when the resulting PeV protons collide with infrared photons of the dust echo (Winter & Lunardini 2023).A potential challenge to a jet-based particle acceleration mechanisms is discussed below in section 5.5.

Contribution to the high-energy neutrino flux
About half of the neutrinos in the IceCube alerts dataset are expected to be background (i.e., atmospheric) events.Based on the "signalness" probability (Aartsen et al. 2017a) of the neutrinos in our sample, we expect about 16 astrophysical neutrinos in our sample.Hence three coincident events implies that at least 19 +22 −12 % (90%CL) of the IceCube astrophysical neutrino alerts are explained by accretion flares in our sample.The fraction of the total astrophysical high-energy neutrino flux produced by the entire accretion flare population is larger than this estimate because our sample of flares is not complete.A factor ≈ 2 increase could be expected to account for neutrino alerts from flares that are too distant to be detected by ZTF and NEOWISE (Stein et al. 2021).
The use of dust echo properties to define the accretion flare population could provide another source of incompleteness.In this work, we use the echoes as tracers of energetic events in the accretion disk of massive black holes.A causal relation between the echo and the neutrino is not required for our analysis, but our search will, by construction, not identify neutrinos from TDEs in dust-free galaxies.

Efficiency puzzle
A sizable contribution of accretion flares to the high-energy neutrino flux is remarkable because the energy we received from these flares is much lower compared to regular AGN.We can estimate the difference in fluence (energy received at Earth) of the two populations using the NEOWISE infrared observations.The sum of the infrared echo flux of the 63 accretion flares is 0.1 Jy, or a fluence of 10 20.5 erg cm −2 for a one year duration of the flare.While the sum of the baseline infrared flux of the AGN detected by ZTF is 26 Jy, a fluence of 10 23.3 erg cm −2 over the 3 year duration of our search (Fig. 1).Here we have only summed the infrared emission of AGN detected in the ZTF alert stream (due to their variability), if we add the contribution of the rest of the population, this estimate of the fluence could increase with another order of magnitude.
Given that steady AGN outshine accretion flares by at least three orders of magnitude we reach a puzzling conclusion: in order to explain the observed PeV-scale neutrino associations, the accretion flares appear to be vastly more efficient at producing PeV neutrinos compared to the majority of AGN.Here we define efficiency as number of high-energy (∼ PeV) particles relative to the electromagnetic energy output (i.e.,   in Eq. 11).Since the contribution of accretion flares and normal AGN to the total neutrino flux have to add to 100% (or less, if we also include other sources of neutrinos), equal efficiency would imply that accretion flares should produce at most 0.1% of the total high-energy neutrino flux.This is clearly not consistent with the lower limit of 10% based on the three accretion flares that are detected as potential neutrino sources.
The fact that accretion flares are a cosmic minority presents a challenge for models of TDEs as neutrino sources that involve a relativistic jet (Farrar & Gruzinov 2009;Wang et al. 2011;Winter & Lunardini 2021) or a corona (Murase et al. 2020).Since AGN also have these features, a similar efficiency of PeV-scale neutrino production can be expected for such models.Here we have formulated the efficiency puzzle in terms of high-energy neutrino production.But unless the optical depth for pion production is much lower for accretion flares compared to AGN, the same conclusions hold for the efficiency of high-energy particle acceleration.
-1500 -1000 -500 0 500 Rest-frame days since peak For each source, the neutrino arrived (dotted vertical lines) a few months after the peak of the optical light curve (red and green symbols for the ZTF  and -band, respectively).The infrared light curve (blue and purple symbols for the WISE W1 and W2, respectively) evolves on longer timescales due to the large distance of the dust sublimation radius (∼ 0.1 pc).From the duration of the dust reverberation signal we infer a peak luminosity near the Eddington limit for all three flares (Table 1).
The low black hole mass of the accretion flares coincident with neutrinos points to a potential solution for this efficiency puzzle.Both TDEs and extreme AGN flares are commonly observed to reach the Eddington limit (e.g., Wevers et al. 2017), while the vast majority of AGN accrete an order of magnitude below the Eddington limit (e.g., Kelly & Shen 2013).Due to photon trapping, a high accretion rate increases the scale height of the accretion disk relative to the geometrically thin disk that operates at Eddington ratios of ≈ 10%.As the Eddington ratio approaches unity, we might expect a state change of the accretion disk (Abramowicz & Fragile 2013).In the context of TDEs, this state change is supported by observations: the X-ray spectra of typical AGN are hard and non-thermal, while TDE and extreme AGN flares have thermal and soft X-ray spectra (Saxton et al. 2020;Frederick et al. 2021;Wevers 2020;Wevers et al. 2021).We can speculate that the accretion state that corresponds to high Eddington ratios enables more efficient particle acceleration.This possibility has been explored by Hayasaki & Yamazaki (2019) for a magnetically arrested disk (MAD).The MAD accretion regime (Narayan et al. 2003) has also been employed by Scepi et al. (2021) to explain the properties of a peculiar AGN flare (1ES 1927+654: Trakhtenbrot et al. 2019).
Because the disk environment will absorb the gamma-ray emission produced in  0 decay through the pair-creation process (Hayasaki & Yamazaki 2019;Murase et al. 2020), a disk-based particle accelerator is predicted to be dark above GeV energies, consistent with the upper limits on gamma-ray emission for the three accretion flares with neutrino counterparts (Table 1).
In summary, the detection of three neutrinos from the energetically-subdominant population of accretion flares can be explained if the high-energy particle acceleration efficiency drastically increases towards the Eddington limit.This scenario might also provides an explanation for neutrino emission from NGC 1068, the most significant hotspot in the IceCube sky map at sub-PeV energies, detected at 4 post-trial (Aartsen et al. 2020b;IceCube Collaboration et al. 2022).NGC 1068 is exceptional because it could be the nearest example of the small subset of AGN accreting near the Eddington limit (Kawaguchi 2003;Lodato & Bertin 2003), similar to the accretion flares presented in this work.The small subset of persistent AGN that accrete close to the Eddington luminosity could provide an important contribution to the potential correlation (detected at 2.6) between persistent AGN and IceCube neutrinos (Abbasi et al. 2022).
and Lucille Packard Foundation.This work was supported by the GROWTH project funded by the National Science Foundation under Grant No 1545949.

Figure 1 .
Figure1.Distribution of the time of maximum light of the ZTF light curve and the detection date of the IceCube alerts.The lack of ZTF events with a peak after mid-2020 happens because the NEOWISE data release includes observations up to the end of 2020 and to be able to measure the properties of the dust echo we require at least two post-peak detections in NEOWISE.Given the 6 month cadence of NEOWISE, the ZTF light curve has to peak prior to mid-2020 to meet this requirement.Because we allow the neutrino to arrive after the peak of the optical flare, all IceCube events up to the end of 2020 are included in the analysis.

Figure 2 .
Figure 2. Parameters inferred from a Gaussian rise plus exponential decay model applied to the ZTF light curves.The flux increase is measured relative to the ZTF reference image.The dashed lines indicates the box that is used to separate accretion flares from regular AGN variability.This requirement selects nuclear supernovae plus all spectroscopically-confirmed ZTF TDEs.The label 'TDE?' indicates accretion flares that occurred in active galaxies (i.e., sources with evidence for accretion prior to the main flare).The three events coincident with a high-energy neutrino are indicated with solid symbols.

Figure 3 .
Figure 3.The dust echo flux and luminosity of nuclear transients in ZTF.We see that the three accretion flares coincident with a high-energy neutrino (filled symbols) are among the strongest dust echoes in the sample of nuclear transients.

Figure 4 .Figure 5 .
Figure 4. Optical spectra of the three accretion flares coincident with a high-energy neutrino.The observed time relative to the peak of the optical emission is indicated.Including data from Frederick et al. (2021); Stein et al. (2021); Reusch et al. (2022).

Figure 6 .
Figure6.Significant dust echoes occur almost exclusively in low-mass black holes.The onset of strong echoes, measured using the infrared flux increase within one year of the optical peak of the flare, coincides with the Hills mass(Hills 1975) for a solar-type star (defined by the requirement that the tidal radius is larger than the black hole horizon).The label 'TDE?' indicates accretion flares that occurred in AGN.The three accretion flares coincident with a high-energy neutrino are indicated with filled symbols.

Figure 9 .
Figure 9. Sky maps of ZTF nuclear transients and IceCube neutrino alerts.In each panel, the cyan boxes indicate the IceCube 90%CL reconstruction areas.The top-left map shows the 63 accretion flares with potential dust echoes that are used in this study.The top-right map shows the density distribution of our "parent sample" of 3142 nuclear flares.The middle-left map shows the result of applying a Gaussian kernel density estimation (KDE) to this sample.The middle-right map shows the resulting sky distribution of Monte Carlo samples that are used for the background hypothesis in our likelihood method.These samples are drawn from the KDE, with the additional exclusion of the Galactic plane (|| > 8, as shown by the grey line).A linear color scale is used in all these panels.The three figures on the bottom row show close ups of the location of the three accretion flares coincident with an IceCube alert.

Figure 10 .
Figure10.The distribution of the test statistics after redistributing the accretion flares in the ZTF sky.From this Monte Carlo simulation we find a probability  = 1.5 × 10 −4 that the three accretion flares we associate with a high-energy neutrino are explained by a chance coincidence.The peak at zero is due to Monte Carlo realizations that yield zero coincident events.

Figure 12 .
Figure12.Cumulative distribution of the dust echo flux.Here we check the assumption that  (Δ IR |) ∝ Δ IR for the signal hypothesis of the likelihood method (Eq.8).For each of the 18 TDE candidates (i.e., accretion flares with  BH < 10 8  ⊙ ) we predict the neutrino expectation value using a linear coupling with echo flux (normalized to yield 1 neutrino for a echo flux of 0.1 mJy).By drawing from a Poisson distribution we simulate neutrino detections from this population.The predicted echo flux distribution for samples that yield three neutrino detections from three different flares (cyan curve) is consistent with the observed flux of the three neutrino-detected flares (black dashed line).

Figure 13 .
Figure13.Neutrino detections for three accretion flares.For each source, the neutrino arrived (dotted vertical lines) a few months after the peak of the optical light curve (red and green symbols for the ZTF  and -band, respectively).The infrared light curve (blue and purple symbols for the WISE W1 and W2, respectively) evolves on longer timescales due to the large distance of the dust sublimation radius (∼ 0.1 pc).From the duration of the dust reverberation signal we infer a peak luminosity near the Eddington limit for all three flares (Table1).

Table 1 .
Multi-messenger inference.The time difference of the neutrino arrival (Δ) is measured relative to the optical peak of the light curve.The angular offset (Δ) is measured relative to the best-fit neutrino position (Table