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

ABSTRACT

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 per cent 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.

neutrinos, galaxies: active, transients: tidal disruption events

1 INTRODUCTION

Accreting black holes have long been suggested as potential sources of high-energy particles (Farrar & Gruzinov 2009; Gaisser & Karle 2017) and this expectation was supported by the detection of a high-energy neutrino coincident (at the 3σ level) with gamma-ray flaring from the blazar TXS 0506+056 (IceCube Collaboration 2018) and from the nearby active galactic nucleus (AGN) NGC 1068 at the 4σ level (IceCube Collaboration 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, 2020a; Hooper, Linden & Vieregg 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 centre 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⁴ AGN detected by ZTF (van Velzen et al. 2021b), less than 1 per cent 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 Near-Earth Object Wide-field Infrared Survey Explorer (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 (T ∼ 2 × 10³ K) at distances of 0.1–1 pc from the black hole (Lu, Kumar & Evans 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 (R_s ∼ 0.1 pc) will yield a similar-looking dust echo: a flat-topped light curve with a duration of 2R_s/c. 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.¹

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.

2 CATALOGUE 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 IceCube neutrino sample and our definition for a flare-neutrino coincidence.

2.1 ZTF nuclear flare selection

As described in van Velzen et al. (2019, 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 host-flare offset (van Velzen et al. 2019) >0.5 arcsec. 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).

Figure 1.

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.

To measure the peak flux of the ZTF light curve, we use the observation with the highest flux, after restricting to 90 per cent of the data points with the highest signal-to-noise ratio (excluding 10 per cent 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-time-scale and fade-time-scale of the ZTF light curve (see Fig. 2).These time-scales are obtained from the measurements of van Velzen et al. (2021b) who applied a Gaussian-rise exponential-decay model to the ZTF alert photometry (both the g and r bands 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) time-scale of the majority of the data points. To remove regular variability from normal AGN, we require a minimal amplitude of the flare of Δm < −1 (m denotes the magnitude) and also we set an upper limit to the rise and fade time-scale (e-folding time <150 and <500 d, 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 (Frederick et al. 2021; Hammerstein et al. 2021; van Velzen et al. 2021b). About 15 per cent of the sources are spectroscopically confirmed SN.

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.

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).

2.2 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 h period of this visit (Wright et al. 2010). The inverse-variance weighted mean of the catalogued flux during these visits is used to construct the NEOWISE 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).

To measure the echo strength we require two observations after the ZTF peak. We define the dust echo flux as ΔF_IR: the difference between the baseline flux and the mean NEOWISE flux within 1 yr after the ZTF peak. The echo strength is ΔF_IR/F_rms, where F_rms is the root-mean-square variability of the baseline observations. The significance of the rms variability is measured using the ratio F_rms/σ_F, where σ_F 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: ΔF_IR/F_rms > F_rms/σ_F. 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 1). In Fig. 3, we show the echo flux and luminosity (for the subset of sources with a spectroscopic redshift) versus echo strength.

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.

Table 1.

Optical and infrared properties of the accretion flares used in this work.

ZTF name	|$t_{\rm peak}^{a}$|	Riseb	Fadeb	ΔFIR/Frms	ΔFIR	z	MBH	|$\frac{P({\rm AGN})}{P(\rm TDE)}^{c}$|	Spectro.
	(MJD)	(d)	(d)		(mJy)		(log10M⊙)		class
AT2019dsg	58620.2	32.1	81.9	92.2	1.58 ± 0.02	0.0512	6.74 (C21)	0.000	TDE
AT2019fdr	58672.5	30.8	336.6	39.2	0.71 ± 0.07	0.2666	7.10 (F20)	0.000	TDE?
AT2019aalc	58658.2	49.0	167.7	15.7	11.13 ± 0.10	0.0356	7.23 (L19)	0.146	TDE?
AT2018dyk	58261.4	60.0	342.0	23.8	1.41 ± 0.03	0.0367	5.50 (F19)	0.000	TDE?
AT2019aame	58363.2	–	138.0	12.3	0.18 ± 0.01	–	–	0.540	–
AT2018lzs	58378.2	134.0	15.5	3.3	0.03 ± 0.01	–	–	3.651	–
AT2021aetz	58390.3	8.6	–	47.5	0.81 ± 0.01	0.0879	6.21 (T13)	0.000	TDE?
AT2018ige	58432.5	–	67.9	65.8	0.31 ± 0.01	–	–	0.000	–
AT2021aeud	58448.3	137.6	282.7	6.2	0.15 ± 0.01	–	–	3.937	–
AT2018iql	58449.4	17.9	41.0	30.1	0.48 ± 0.01	–	–	0.000	–
AT2018jut	58449.6	–	–	5.0	0.12 ± 0.01	–	–	5.424	–
AT2021aeue	58475.1	48.5	48.8	4.9	0.11 ± 0.01	–	–	5.552	–
AT2019aamf	58506.4	80.0	141.2	6.6	0.19 ± 0.01	–	–	3.395	–
AT2018kox	58510.2	26.9	167.6	5.6	0.33 ± 0.01	0.096	–	4.668	TDE?
AT2018lhv	58513.5	17.5	116.1	32.3	0.35 ± 0.01	–	–	0.000	–
AT2019avd	58534.3	14.9	52.9	67.5	1.87 ± 0.05	0.0296	6.10 (F20)	0.000	TDE?
AT2016eix	58539.4	104.2	43.1	6.9	0.16 ± 0.01	–	–	2.922	–
AT2019aamg	58540.5	–	93.7	8.3	0.14 ± 0.00	–	–	1.612	–
AT2018lcp	58547.2	95.7	297.4	12.7	0.16 ± 0.01	0.06	–	0.482	TDE?
AT2021aeuf	58556.4	17.0	137.5	15.6	0.18 ± 0.01	–	–	0.155	–
AT2020aezy	58558.4	82.5	359.8	4.8	0.06 ± 0.01	–	–	5.623	–
AT2019cle	58568.4	9.2	54.8	19.4	0.25 ± 0.01	–	–	0.012	–
AT2019aamh	58582.5	–	154.5	7.7	0.40 ± 0.01	–	–	2.095	–
AT2019dll	58605.2	29.6	–	6.8	0.26 ± 0.01	0.101	7.48 (T13)	3.109	TDE?
AT2019gur	58607.5	–	234.5	38.6	0.34 ± 0.01	–	–	0.000	–
AT2018lof	58608.2	70.9	370.4	4.1	0.15 ± 0.01	0.302	8.98 (L19)	6.078	AGN
AT2019dqv	58628.2	67.0	475.0	40.4	1.77 ± 0.02	0.0816	6.67 (L19)	0.000	TDE?
AT2019cyq	58637.2	49.5	95.0	31.8	0.54 ± 0.01	0.262	7.56 (L19)	0.000	TDE?
AT2021aeug	58641.2	48.6	68.8	4.6	0.11 ± 0.01	–	–	5.900	–
AT2019ihv	58646.5	–	19.3	8.7	0.29 ± 0.01	0.1602	–	1.395	–
AT2019dzh	58651.2	51.4	348.9	6.4	0.15 ± 0.01	0.314	–	3.659	–
AT2019kqu	58652.2	95.6	315.6	6.1	0.27 ± 0.01	0.174	7.53 (L19)	3.996	TDE?
AT2019hbh	58652.3	19.8	106.7	8.4	0.75 ± 0.01	–	–	1.556	–
AT2020aezz	58677.3	99.8	280.6	5.8	0.18 ± 0.01	–	–	4.405	–
AT2020afaa	58678.2	83.8	330.3	7.0	0.13 ± 0.02	–	–	2.830	–
AT2019idm	58682.2	51.7	–	25.2	0.41 ± 0.02	0.0544	6.64 (T13)	0.000	TDE?
AT2019ihu	58709.5	79.7	110.9	6.2	0.25 ± 0.02	0.27	8.90 (L19)	3.843	AGN
AT2019meh	58713.1	23.2	127.9	29.7	1.99 ± 0.04	0.0935	7.06 (V23)	0.000	TDE?
AT2020afab	58717.2	47.5	127.2	5.0	0.22 ± 0.01	0.2875	6.49 (V23)	5.432	TDE?
AT2019aami	58717.4	–	152.0	31.8	0.41 ± 0.01	–	–	0.000	–
AT2019nna	58717.4	34.4	87.1	27.0	0.36 ± 0.01	–	–	0.000	–
AT2019nni	58732.2	29.1	109.5	4.9	0.36 ± 0.01	0.137	–	5.489	–
AT2021aeuk	58733.1	43.4	238.1	7.3	0.33 ± 0.01	0.235	–	2.464	–
AT2019hdy	58749.5	–	87.9	4.0	0.20 ± 0.01	0.442	–	6.000	–
AT2019pev	58750.1	13.1	54.7	7.4	0.19 ± 0.01	0.097	6.40 (F20)	2.364	TDE?
AT2019qiz	58753.1	11.9	34.7	44.4	2.14 ± 0.04	0.01499	6.19 (N20)	0.000	TDE
AT2019brs	58758.1	110.9	483.2	9.6	0.78 ± 0.01	0.3736	7.20 (F20)	1.060	TDE?
AT2020afac	58758.3	68.0	95.0	10.8	0.16 ± 0.01	–	–	0.802	–
AT2019wrd	58764.3	82.6	–	7.6	0.15 ± 0.01	–	–	2.197	–
AT2021aeuh	58789.5	–	–	3.9	0.28 ± 0.01	0.0834	7.39 (L19)	5.867	TDE?
AT2019msq	58791.2	137.9	–	6.4	0.16 ± 0.01	–	–	3.586	–
AT2019qpt	58798.3	44.5	135.1	13.8	0.24 ± 0.01	0.242	6.97 (L19)	0.341	TDE?
AT2020afad	58802.2	77.1	125.5	3.7	0.08 ± 0.01	–	–	5.216	–
AT2019mss	58811.6	35.2	209.7	20.8	0.38 ± 0.01	–	–	0.004	–
AT2019thh	58851.1	104.7	433.8	72.2	1.75 ± 0.02	0.0506	–	0.000	TDE?
AT2021aeui	58860.3	–	60.0	6.2	0.49 ± 0.02	–	–	3.899	–
AT2020afae	58867.2	30.2	42.0	5.2	0.35 ± 0.02	–	–	5.200	–
AT2020mw	58867.3	12.6	32.1	6.8	0.12 ± 0.01	–	–	3.061	–
AT2020iq	58878.1	22.7	71.5	24.5	0.45 ± 0.01	0.096	6.37 (V23)	0.000	TDE?
AT2019xgg	58891.2	81.6	446.2	4.4	0.13 ± 0.01	–	–	6.003	–
AT2020atq	58903.2	39.3	205.5	20.8	0.66 ± 0.01	–	–	0.003	–
AT2021aeuj	58974.2	79.0	233.5	18.1	0.20 ± 0.01	0.695	–	0.032	–
AT2020hle	58978.3	32.7	62.9	21.0	0.33 ± 0.01	0.103	6.40 (F20)	0.003	TDE?

Note. The first three entries of this table list the events coincident with an IceCube neutrino alert. ^aThe column t_peak lists the time of maximum light of the ZTF light curve. ^bThe rise and fade columns list the e-folding time. ^cThe ratio of the AGN and TDE probability is based on the strength of the dust echo (ΔF_IR/F_rms, see Fig. 7). The references listed behind the black hole mass estimates give the origin of the optical spectrum that was used: T13 (Thomas et al. 2013), L19 (Liu et al. 2019), F19 (Frederick et al. 2019), F20 (Frederick et al. 2021), N20 (Nicholl et al. 2020), C21 (Cannizzaro et al. 2021), V23 (this work).

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, Kumar & Evans 2016; van Velzen et al. 2016, 2021a). Using this geometric luminosity estimate (equation 12 in van Velzen et al. 2016), we find a bolometric luminosity |$\sim 10^{45}\, {\rm erg}\, {\rm s}^{-1}$| for all accretion flares coincident with high-energy neutrinos. All are consistent with reaching the Eddington limit (Table 3). 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.

Table 3.

Multimessenger inference. The time difference of the neutrino arrival (Δt) is measured relative to the optical peak of the light curve. The angular offset (Δd) is measured relative to the best-fitting neutrino position (Table 2). The upper limit on the neutrino luminosity is estimated by assuming an expectation value of one particle in Δt (the true neutrino luminosity could be one or two orders of magnitude lower due to the Eddington bias, see Section 5.2). The X-ray temperature (T_X) and luminosity (L_X) are based on SRG/eROSITA or Swift/XRT and cover an energy range of 0.2 to 10 keV (see Section 4.3). The radio luminosity is measured at ν ≈ 5 GHz (see Section 4.1). The 95 per cent CL upper limit on the gamma-ray luminosity (L_γ) is obtained from data of the Fermi-LAT telescope and covers photons in the energy range 0.1–800 GeV (see Section 4.2). The black hole mass (M_BH) is estimated from the optical spectrum of the host galaxy (Section 2.3). The Eddington ratio (f_Edd) follows from the bolometric luminosity (L_bol) as estimated from the duration of the dust echo (Section 2.2).

Flare	Neutrino	z	Δt	Δd	Lneutrino	TX	LX	Lradio	Lγ	Lbol	MBH	fEdd
			(d)	(deg)	(⁠|${\rm erg}\, {\rm s}^{-1}$|⁠)	(keV)	(⁠|${\rm erg}\, {\rm s}^{-1}$|⁠)	(⁠|${\rm erg}\, {\rm s}^{-1}$|⁠)	|$({\rm erg}\, {\rm s}^{-1}$|⁠)	(⁠|${\rm erg}\, {\rm s}^{-1}$|⁠)	(M⊙)
AT2019dsg	IC191001A	0.051	150	1.3	<1044.5	0.07 ± 0.01	1043.4	1038.9	<1043.1	1045.3	106.6	4.3
AT2019fdr	IC200530A	0.267	289	1.7	<1045.8	0.06 ± 0.03	1043.1	1039.3	<1044.3	1044.9	107.1	0.5
AT2019aalc	IC191119A	0.036	148	1.9	<1044.2	0.17 ± 0.01	1042.1	1038.7	<1043.3	1045.1	107.2	0.6

2.3 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 M–σ 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 M–σ 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 disc luminosity. We adopt the relation from Ho & Kim (2015):

$$\begin{eqnarray} \log (M_{\rm BH} / {\rm M}_\odot) &=& \log \Bigg [\left(\frac{{\rm FWHM~(H}\,\beta)}{\rm 1000~km~s^{-1}}\right)^2 \\ &&\times \left(\frac{L_{5100}}{10^{44}\, {\rm erg}\, {\rm s}^{-1}}\right)^{0.533} \Bigg ] + 6.91. \end{eqnarray}$$

(1)

Here, L₅₁₀₀ 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 catalogue of Liu et al. (2019), who selected 14 584 type-1 AGN based on detection of a broad H α line and applied equation (1) to measure the black hole mass. We find 580 nuclear transients with black hole mass estimates based on this catalogue. 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, 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 2 yr after the peak of the optical flare.

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), and Reusch et al. (2022).

The spectra of AT2020iq, AT2019meh, and AT2019aalc were obtained with the Low Resolution Imaging Spectrograph (LRIS; Oke et al. 1995) on the 10-m Keck-I telescope at 20, 660, and 714 d post t_peak (Table 1), respectively. The new spectrum of AT2020afab was obtained with the Double Spectrograph (DBSP; Oke & Gunn 1982) on the 5-m Palomar telescope (P200) at 15 d post t_peak. The DBSP spectrum was reduced with the pyraf-dbsp pipeline (Bellm & Sesar 2016). The LRIS spectra were reduced using lpipe (Perley 2019).

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) M–σ relation. This yields 219 additional M_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 1, we list the reference for the black hole mass estimate for each accretion flare.

2.4 IceCube neutrino alerts

Our parent sample of neutrino alerts (Aartsen et al. 2017a) includes the events published up to 2020 December 31. We exclude the IceCube 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 per centCL area larger than 300 deg² (IC200410A; Stein 2020c). This leaves 43 events, listed in Table 2. Of these, three fall outside ZTF extra-galactic footprint (Ω_ZTF = 2.8 × 10⁴ deg²; Stein et al. 2021), leaving 40 IceCube alerts that can yield a coincidence with an accretion flare in our ZTF+NEOWISE data set.

Table 2.

IceCube neutrino alerts used in this study.

Event	RA	Decl.	90 per cent CL area	Signalness	GCN ref.
			(deg2)
IC180908A	|$144.58_{-1.55}^{+1.45}$|	|$-2.13_{-0.90}^{+1.20}$|	6.3	0.34	Blaufuss (2018a)
IC181023A	|$270.18_{-2.00}^{+1.70}$|	|$-8.57_{-1.25}^{+1.30}$|	9.3	0.28	Blaufuss (2018b)
IC190704A	|$161.85_{-4.33}^{+2.16}$|	|$27.11_{-1.83}^{+1.81}$|	21.0	0.49	Santander (2019a)
IC190712A	|$76.46_{-6.83}^{+5.09}$|	|$13.06_{-3.44}^{+4.48}$|	92.1	0.30	Blaufuss (2019c)
IC190819A	|$148.80_{-3.24}^{+2.07}$|	|$1.38_{-0.75}^{+1.00}$|	9.3	0.29	Santander (2019b)
IC190922A	|$167.43_{-2.63}^{+3.40}$|	|$-22.39_{-2.89}^{+2.88}$|	32.2	0.51	Stein (2019b)
IC191119A	|$230.10_{-6.48}^{+4.76}$|	|$3.17_{-2.09}^{+3.36}$|	61.2	0.45	Blaufuss (2019e)
IC191122A	|$27.25_{-2.90}^{+1.70}$|	|$-0.04_{-1.49}^{+1.17}$|	12.2	0.33	Blaufuss (2019f)
IC191204A	|$79.72_{-1.74}^{+3.20}$|	|$2.80_{-1.23}^{+1.12}$|	11.6	0.33	Stein (2019d)
IC191215A	|$285.87_{-3.19}^{+2.88}$|	|$58.92_{-2.25}^{+1.85}$|	12.8	0.47	Stein (2019e)
IC191231A	|$46.36_{-3.47}^{+4.27}$|	|$20.42_{-2.80}^{+2.11}$|	35.5	0.46	Santander (2019c)
IC200421A	|$87.93_{-2.83}^{+3.44}$|	|$8.23_{-1.84}^{+2.09}$|	24.4	0.33	Blaufuss (2020a)
IC200425A	|$100.10_{-3.14}^{+4.67}$|	|$53.57_{-1.60}^{+2.45}$|	19.0	0.48	Santander (2020a)
IC200523A	|$338.64_{-6.07}^{+10.77}$|	|$1.75_{-3.54}^{+1.84}$|	90.5	0.25	Blaufuss (2020b)
IC200614A	|$33.84_{-6.39}^{+4.77}$|	|$31.61_{-2.28}^{+2.75}$|	47.9	0.42	Blaufuss (2020c)
IC200615A	|$142.95_{-1.45}^{+1.18}$|	|$3.66_{-1.06}^{+1.19}$|	5.9	0.83	Lagunas Gualda (2020c)
IC200806A	|$157.25_{-0.89}^{+1.21}$|	|$47.75_{-0.64}^{+0.65}$|	1.8	0.40	Stein (2020e)
IC200911A	|$51.11_{-11.01}^{+4.42}$|	|$38.11_{-1.99}^{+2.35}$|	52.8	0.41	Lagunas Gualda (2020d)
IC200921A	|$195.29_{-1.73}^{+2.35}$|	|$26.24_{-1.77}^{+1.51}$|	12.0	0.41	Lagunas Gualda (2020e)
IC200926B	|$184.75_{-1.55}^{+3.64}$|	|$32.93_{-0.91}^{+1.15}$|	9.0	0.43	Blaufuss (2020e)
IC201014A	|$221.22_{-0.75}^{+1.00}$|	|$14.44_{-0.46}^{+0.67}$|	1.9	0.41	Blaufuss (2020f)
IC201114A	|$105.25_{-1.12}^{+1.28}$|	|$6.05_{-0.95}^{+0.95}$|	4.5	0.56	Blaufuss (2020g)
IC201115A	|$195.12_{-1.49}^{+1.27}$|	|$1.38_{-1.11}^{+1.30}$|	6.6	0.49	Lagunas Gualda (2020i)
IC201120A	|$307.53_{-5.59}^{+5.34}$|	|$40.77_{-2.80}^{+4.97}$|	65.3	0.50	Lagunas Gualda (2020j)
IC201221A	|$261.69_{-2.50}^{+2.29}$|	|$41.81_{-1.20}^{+1.29}$|	8.9	0.56	Blaufuss (2020h)
IC190503A	|$120.28_{-0.77}^{+0.57}$|	|$6.35_{-0.70}^{+0.76}$|	1.9	0.36	Blaufuss (2019a)
IC190619A	|$343.26_{-2.63}^{+4.08}$|	|$10.73_{-2.61}^{+1.51}$|	27.1	0.55	Blaufuss (2019b)
IC190730A	|$225.79_{-1.43}^{+1.28}$|	|$10.47_{-0.89}^{+1.14}$|	5.4	0.67	Stein (2019a)
IC190922B	|$5.76_{-1.37}^{+1.19}$|	|$-1.57_{-0.82}^{+0.93}$|	4.5	0.51	Blaufuss (2019d)
IC191001A	|$314.08_{-2.26}^{+6.56}$|	|$12.94_{-1.47}^{+1.50}$|	25.5	0.59	Stein (2019c)
IC200109A	|$164.49_{-4.19}^{+4.94}$|	|$11.87_{-1.36}^{+1.16}$|	22.5	0.77	Stein (2020b)
IC200117A	|$116.24_{-1.24}^{+0.71}$|	|$29.14_{-0.78}^{+0.90}$|	2.9	0.38	Lagunas Gualda (2020a)
IC200512A	|$295.18_{-2.26}^{+1.72}$|	|$15.79_{-1.29}^{+1.26}$|	9.8	0.32	Lagunas Gualda (2020b)
IC200530A	|$255.37_{-2.56}^{+2.48}$|	|$26.61_{-3.28}^{+2.33}$|	25.2	0.59	Stein (2020d)
IC200620A	|$162.11_{-0.95}^{+0.64}$|	|$11.95_{-0.48}^{+0.63}$|	1.7	0.32	Santander (2020b)
IC200916A	|$109.78_{-1.44}^{+1.08}$|	|$14.36_{-0.85}^{+0.88}$|	4.2	0.32	Blaufuss (2020d)
IC200926A	|$96.40_{-0.55}^{+0.73}$|	|$-4.33_{-0.76}^{+0.61}$|	1.7	0.43	Lagunas Gualda (2020f)
IC200929A	|$29.53_{-0.53}^{+0.53}$|	|$3.47_{-0.35}^{+0.71}$|	1.1	0.47	Lagunas Gualda (2020g)
IC201007A	|$265.17_{-0.52}^{+0.52}$|	|$5.34_{-0.23}^{+0.32}$|	0.6	0.89	Santander (2020c)
IC201021A	|$260.82_{-1.68}^{+1.73}$|	|$14.55_{-0.74}^{+1.35}$|	6.9	0.30	Lagunas Gualda (2020h)
IC201130A	|$30.54_{-1.31}^{+1.13}$|	|$-12.10_{-1.13}^{+1.15}$|	5.4	0.15	Lagunas Gualda (2020k)
IC201209A	|$6.86_{-1.22}^{+1.02}$|	|$-9.25_{-1.14}^{+0.99}$|	4.7	0.19	Lagunas Gualda (2020l)
IC201222A	|$206.37_{-0.80}^{+0.90}$|	|$13.44_{-0.38}^{+0.55}$|	1.5	0.53	Blaufuss (2020i)

The signalness, or P_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 per cent of neutrinos whose true location is outside the reported 90 per centCL area, we estimate that about 16 cosmic neutrinos can in principle be recovered by our analysis.

2.5 Flare-neutrino coincidence

An accretion flare is considered coincident with a neutrino when the source falls inside the 90 per centCL 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 data set 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 echoes compared to the rest of the nuclear flare population (Fig. 5). The significance of this result is discussed below in Section 3.

Figure 5.

Flares coincident with high-energy neutrinos have very strong dust echoes. The TDE candidates (blue line) are defined as accretion flares with an estimated black hole mass (Fig. 6) below |$10^{8}\, M_\odot$|⁠. This population shows much stronger dust echoes compared to regular variability from the AGN population (grey line). In the full population of 140 nuclear transients with candidate dust echoes (pink line), we find three events coincident with a high-energy neutrino. The echo strength of these three flares is inconsistent with originating from the full population (p = 0.007 for an Anderson–Darling test).

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 (ΔF_IR/F_rms > 10) are observed almost exclusively for events with |$M_{\rm BH}\lt 10^{8}\, {\rm M}_\odot$| (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 |$M_{\rm BH}\lt 10^{8}\, {\rm M}_\odot$|⁠. 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.

Figure 6.

Significant dust echoes occur almost exclusively in low-mass black holes. The onset of strong echoes, measured using the infrared flux increase within 1 yr 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.

3 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 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).

3.1 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 several possible hypotheses that can be used to explain the origin. If there are N accretion flares, we can test N discrete hypotheses, i.e. that the neutrino originated from source #1 (H₁), from source #2 (H₂), etc. We can denote this group of hypotheses as signal hypotheses, i.e. H_i ≡ S_i. 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 H₀. We can also denote this null hypothesis, or background hypothesis, as H₀ ≡ B.

For each hypothesis, we can define the likelihood that the data, x, is obtained under that hypothesis |$\mathcal {L}_{i} = P(x | H_{i})$|⁠. For each neutrino, we select the hypothesis with the greatest likelihood as our best fit:

$$\begin{eqnarray} \hat{\mathcal {L}} = \max \left(\mathcal {L}_{i} \right). \end{eqnarray}$$

(2)

In the rare occasion that multiple accretion flares are found in coincidence with a single neutrino, we select the flare with the highest likelihood:

$$\begin{eqnarray} \hat{\mathcal {L}_{s}} = {\rm max} \left(\mathcal {L}_{i\gt 0} \right). \end{eqnarray}$$

(3)

Next, we compare the signal hypothesis to the background hypothesis:

$$\begin{eqnarray} \hat{\mathcal {L}} = {\rm max} \left(\hat{\mathcal {L}_{s}}, \mathcal {L}_{0} \right). \end{eqnarray}$$

(4)

When no accretion flare is found in coincidence with a neutrino, we have |$\hat{\mathcal {L}} =\mathcal {L}_0$|⁠.

Below, we first discuss the input for the signal hypothesis, followed by the components of the background hypothesis.

3.1.1 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:

• 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 P_coin(α, δ, t|S), with α and δ the right ascension and declination of the accretion flare, respectively.

• 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 per cent CL localization area of the neutrino. Hence, for coincident events, P_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 |$\lt 10^8\, M_\odot$|⁠). 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 use of echo flux in the signal likelihood is motivated by applying what is arguably the simplest-possible model for neutrino 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:

$$\begin{eqnarray} P(\Delta F_{\rm IR}|S) = \frac{\Delta F_{\rm IR}}{{\rm max}(\Delta F_{\rm IR})-{\rm min}(\Delta F_{\rm IR})}. \end{eqnarray}$$

(5)

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

$$\begin{eqnarray} P(x | S) = 0.9 \times P(\Delta F_{\rm IR}/F_{\rm rms}|S) \times P(\Delta F_{\rm IR} | S) . \end{eqnarray}$$

(6)

3.1.2 Components of the background hypothesis

We now consider the probability of observing the data under the null hypothesis. We again have two distinct factors:

• The probability of a spatial and temporal coincidence if the neutrino was not produced by any accretion flare: P_coin(α, δ, t|B).

• 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 coincidences, both spatially and temporally. The expectation value for the number of chance coincidences within the 90 per cent CL area (Ω_ν) of a given neutrino is given by n_bg × Ω_ν, with n_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

$$\begin{eqnarray} P_{\rm coin}(\alpha ,\delta , t|B)=n_{\rm bg} ~ \Omega _\nu . \end{eqnarray}$$

(7)

A Monte Carlo approach is used to estimate n_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 transients. 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 |b| < 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⁶ samples of 63 accretion flares, we find a total of 4.7 × 10⁵ 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², the total neutrino area that overlaps the ZTF footprint: n_bg = 0.48/ΣΩ_ν = 6.9 × 10⁻⁴ deg⁻².

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 P(ΔF_IR/F_rms|B). The resulting PDFs are displayed in Figs 7 and 8.

Figure 7.

The observed distribution of dust echo strength for the full population of nuclear transients with variable infrared emission (green). The observed distributions are approximated using kernel density estimation (dotted lines). The ratio of these two probability density functions (solid black line) is used in the likelihood analysis (equation 8).

Figure 8.

The observed distribution of echo flux for the population of nuclear transients with post-peak NEOWISE detections (green). From this parent sample, a PDF is obtained using kernel density estimation (dotted lines). For the signal hypothesis we assume |$P(\rm \Delta F_{\rm IR}|S)\propto \Delta F_{\rm 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).

3.1.3 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 (P_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:

$$\begin{eqnarray} {\rm TS}\, &=& -2\ln \Bigg (\frac{ \hat{\mathcal {L}_{s}}}{\mathcal {L}_{0}}\Bigg) \\ &=&-2 \ln \Bigg (\prod _i P_{\rm astro} \times \frac{0.9}{\Omega _\nu n_{\rm bg}} \\ &~&\times \frac{P(\Delta F_{\rm IR}/F_{\rm rms} | S)}{P(\Delta F_{\rm IR}/F_{\rm rms} | B)} \times \frac{P(\Delta F_{\rm IR} | S)}{P(\Delta F_{\rm IR}|B)} \Bigg) . \end{eqnarray}$$

(8)

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⁶ 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⁻⁴ of the simulations for the background hypothesis have an equal or greater TS than the observations, corresponding to a significance of 3.6σ.

Figure 9.

Sky maps of ZTF nuclear transients and IceCube neutrino alerts. In each panel, the cyan boxes indicate the IceCube 90 per centCL 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 (|b| > 8, as shown by the grey line). A linear colour 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.

The distribution of the test statistics after redistributing the accretion flares in the ZTF sky. From this Monte Carlo simulation we find a probability p = 1.5 × 10⁻⁴ 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.

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 of 60; if our search had not uncovered this new event, the significance would have been 2.4σ.

3.2 Statistical power

It can be instructive to investigate the statistical power of our likelihood method (equation 8) as a function of the expectation value of the number of signal neutrinos, n_exp. To simulate a signal for a given n_exp, we draw M events from a Poisson distribution with expectation value n_exp. Each of these M events has a 90 per cent 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 and 8). We repeat this signal simulation 10 000 times to obtain the TS distribution for a given n_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 per cent 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 n_exp > 0.3. For reference, our detection of three coincident neutrinos implies |$n_{\rm exp}=3_{-1.9}^{+3.8}$| (90 per cent CL). From Fig. 11, we see that for this range of n_exp, a significance between 3 and 4σ should be expected.

Figure 11.

Assessment of the statistical power of our weighted likelihood ratio test (equation 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 (n_exp in equation 9). For each expectation value, one can simulate a distribution of the test statistic (TS; equation 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 catalogue (based on equation 9). For reference, based on the detection of three coincident neutrino-flare pairs, the 90 per centCL range of n_exp is (1.1, 6.8).

The number of signal neutrinos is related to the fraction of cosmic high-energy neutrinos that originate from accretion flares, f_flare:

$$\begin{eqnarray} f_{\rm ac} = \frac{n_{\rm exp}}{N_{\rm cosmic}\eta _{\rm ZTF}}. \end{eqnarray}$$

(9)

Here, N_cosmic ≈ 16 is the number of astrophysical neutrino alerts 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 catalogue. In particular, a large population of relatively high-redshift accretion flares (z > 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 f_ac is shown at the top of Fig. 11.

3.3 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 (n_ν) and the infrared flux: n_ν = ϵΔF_IR, with ϵ = 100 Jy⁻¹. Applying this coupling to the observed echo flux distribution of the TDE candidates, we obtain n_ν 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 per cent 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 per cent 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.

Figure 12.

Cumulative distribution of the dust echo flux. Here, we check the assumption that P(ΔF_IR|S) ∝ ΔF_IR for the signal hypothesis of the likelihood method (equation 8). For each of the 18 TDE candidates (i.e. accretion flares with |$M_{\rm BH}\lt 10^8\, M_\odot$|⁠), we predict the neutrino expectation value using a linear coupling with echo flux (normalized to yield 1 neutrino for an 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).

While this agreement is encouraging, we should expect deviations from some of the simplifying assumptions that are made to construct the test statistic (equation 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 equation (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 P(ΔF_IR|S). This approach is not ideal for two reasons: (1) the coupling strength is a free parameter and (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 p = 1.7 × 10⁻⁴ 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⁴ 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 n_bg = 8.8 × 10⁻⁴ and p = 2.1 × 10⁻⁴ (3.5σ).

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. ΔF_IR/F_rms > F_rms/σ_F), 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_i = 0). The probability of the null hypothesis is p = 3.2 × 10⁻⁴ (3.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 ΔF_IR/F_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 f_temp = 0.24. We thus obtain the effective source density of large echoes: n_eff = 29 × f_temp/Ω_ZTF = 2.5 × 10⁻⁴ deg⁻², with Ω_ZTF = 2.8 × 10⁴ deg² 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⁻⁴ (3.2σ).

4 MULTIWAVELENGTH PROPERTIES

4.1 Detection at radio wavelengths

The first neutrino-detected source AT2019dsg, showed transient radio emission, evolving on a time-scale of month, with a peak luminosity (νL_ν) of |$8 \times 10^{38}\, {\rm erg}\, {\rm s}^{-1}$| at 10 GHz (Stein et al. 2021). The source AT2019fdr was detected in radio follow-up observations with similar luminosity (⁠|$2 \times 10^{39}\, {\rm erg}\, {\rm 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 (Faint Images of the Radio Sky at Twenty; Becker, White & Helfand 1995) or VLASS (Very Large Array Sky Survey; 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\times 10^{38}\, {\rm erg}\, {\rm 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³⁸ and 5 × 10³⁸  |${\rm erg}\, {\rm 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 yr post-peak) is statistically significant (at the 8σ level, as estimated from the rms in the VLASS ‘Quick Look’ images).

4.2 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² region centred 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.fits and iso_P8R3_SOURCE_V2_v1.txt for 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 d 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 catalogue (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 per cent CL upper limit for the energy flux (i.e. integrated over the whole analysis energy range) listed in Table 3 is derived for a power-law spectrum (dN/dE ∝ E^−Γ) with photon power-law index Γ = 2. The Fermi-LAT upper limit for AT2019dsg listed in Table 3 is obtained using the same LAT data analysis strategy and covers a similar time window (150 d) 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).

4.3 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 2021 October, AT2019aalc (= SRGe J152416.7+045118) has been visited four times starting from 2020 February 2, 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⁻¹³ erg s⁻¹cm⁻². The source had a soft thermal spectrum with the best fit blackbody temperature of kT = 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.0 keV flux of ≈6.0 × 10⁻¹⁴ erg s⁻¹cm⁻² and an extremely soft thermal spectrum with a temperature of |$56_{-26}^{+32}$| 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⁻¹⁴ erg s⁻¹cm⁻² (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 3σ upper limit for the combined data of the three visits of ≲ 1.9 × 10⁻¹⁴ erg s⁻¹cm⁻² (assuming a power law spectrum with the slope of Γ = 1.8). The X-ray measurements of AT2019dsg listed in Table 3 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.

5 RESULTS AND IMPLICATIONS

5.1 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 high-energy 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 multiwavelength properties of the three sources with a neutrino counterpart.

For AT2019aalc, we measure a soft thermal spectrum with temperature of kT = 172 ± 10 eV. Such soft thermal emission is rare: of all accretion flares with potential dust echoes, less than 13 per cent are as soft as AT2019aalc (kT < 172 eV), and even fewer are as soft as AT2019dsg and AT2019fdr. These distinctive X-ray properties provide 3σ level (p = 0.13³) 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 per cent 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⁻³.

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 on to massive black holes.

Figure 13.

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 g and r bands, 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 3).

5.2 Energetics

Around PeV energies, the expectation value for the number of IceCube neutrino detections can be approximated as

$$\begin{eqnarray} N_{\rm neutrino} = 0.7 \times \frac{E_{\rm neutrino}}{10^{51}~{\rm erg}} \left(\frac{d}{100~{\rm Mpc}} \right)^{-2}, \end{eqnarray}$$

(10)

with E_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 N_neutrino < 1 for any individual source in our sample. This implies equation (10) cannot be inverted to estimate the energy in neutrinos emitted by a single source. Only a upper limit can be estimated by setting N_neutrino = 1, the result is shown in Table 3.

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 N_neutrino for the entire population, we again assume a single coupling strength (η_ν) between the electromagnetic energy and the energy carried by neutrinos:

$$\begin{eqnarray} E_{\rm neutrino} = \eta _{\nu } E_{\rm bol} \quad . \end{eqnarray}$$

(11)

We estimate the bolometric energy of the flare from its bolometric luminosity as obtained from the duration of the dust echo (L_bol, see Section 2.2) and the e-folding time τ of the decaying part of the optical light curve (see Fig. 2): E_bol = L_bolτ. The total bolometric flare energy of the 63 accretion flares is 1 × 10⁵⁴ erg. Hence, an expectation of |$N_{\rm neutrino} = 3_{-2}^{+4}$| (90 per cent CL) yields a mean coupling strength |$\eta _\nu = 4_{-3}^{+5}\times 10^{-2}$|⁠.

For most models of particle acceleration in AGN accretion discs, the luminosity in high-energy cosmic rays is between one and two orders of magnitude larger than the neutrino luminosity (Begelman, Rudak & Sikora 1990; Murase et al. 2020). Our estimate of the fraction of the electromagnetic energy that is given to neutrinos (η_ν ∼ 10⁻²) 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 disc with a constant accretion rate yields an expectation of 0.1 IceCube high-energy neutrinos for an integration time of 1 yr (Murase et al. 2020). We refer to Winter & Lunardini (2023) for further exploration of potential particle acceleration mechanisms in all three accretion flares.

5.3 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 p = 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 1 yr, 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).

Secondly, the mass accretion rate may not be constant, but delayed. The circularization time-scale of the stellar debris can be estimated as

$$\begin{eqnarray} t_{\rm circ}\approx && 0.55\, \beta ^{-3} \left(\frac{\eta }{0.1}\right)^{-1} \\ &&\left(\frac{M_{\rm bh}}{10^7\, {\rm M}_\odot }\right)^{-7/6} \, \left(\frac{m_*}{1\, {\rm M}_\odot }\right)^{8/25}\, \, \, {\rm yr} \end{eqnarray}$$

(12)

(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, Rossi & Lodato (2017). Because the inner accretion disc is formed after the circularization of the debris, it could take of order t_circ for the first neutrinos to be produced, which appears to match the observed time delay (Fig. 13).

A caveat to equation (12) is a potential interaction of stellar debris stream with a pre-existing accretion disc, which can shorten the circularization time-scale (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 disc, implying the disc has already formed when the optical emission peaks (Bonnerot & Lu 2020). In this case, for PeV particle acceleration in the newly formed accretion disc, 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.

5.4 Contribution to the high-energy neutrino flux

About half of the neutrinos in the IceCube alerts data set 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_{-12}^{+22}{{\ \rm per\ cent}}$| (90 per cent 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 disc 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.

5.5 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}\, {\rm erg}\, {\rm cm}^{-2}$| for a 1 yr 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} \, {\rm erg}\, {\rm cm}^{-2}$| over the 3 yr 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 equation 11). Since the contribution of accretion flares and normal AGN to the total neutrino flux have to add to 100 per cent (or less, if we also include other sources of neutrinos), equal efficiency would imply that accretion flares should produce at most |$0.1{{\ \rm per\ cent}}$| of the total high-energy neutrino flux. This is clearly not consistent with the lower limit of 10 per cent 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.

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 disc relative to the geometrically thin disc that operates at Eddington ratios of ≈10 per cent. As the Eddington ratio approaches unity, we might expect a state change of the accretion disc (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; Wevers 2020; Frederick et al. 2021; 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 disc (MAD). The MAD accretion regime (Narayan, Igumenshchev & Abramowicz 2003) has also been employed by Scepi, Begelman & Dexter (2021) to explain the properties of a peculiar AGN flare (1ES 1927+654: Trakhtenbrot et al. 2019).

Because the disc environment will absorb the gamma-ray emission produced in π⁰ decay through the pair-creation process (Hayasaki & Yamazaki 2019; Murase et al. 2020), a disc-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 3).

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 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).

ACKNOWLEDGEMENTS

We acknowledge useful discussions and suggestions from J. Becerra González, M. Kerr, W. Lu, C. Lunardini, K. Murase, and W. Winter. We thank the anonymous referee for the useful suggestions and comments.

Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the ZTF project. ZTF is supported by the National Science Foundation under grant no. AST-1440341 and grant no. AST-2034437 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, and IN2P3, France. Operations are conducted by COO, IPAC, and UW. SED Machine is based upon work supported by the National Science Foundation under grant no. 1106171.

This work is based on observations with the eROSITA telescope on board the SRG observatory. The SRG observatory was built by Roskosmos in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI) in the framework of the Russian Federal Space Program, with the participation of the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG/eROSITA X-ray telescope was built by a consortium of German Institutes led by MPE, and supported by DLR. The SRG spacecraft was designed, built, launched, and is operated by the Lavochkin Association and its subcontractors. The science data are downlinked via the Deep Space Network Antennae in Bear Lakes, Ussurijsk, and Baykonur, funded by Roskosmos. The eROSITA data used in this work were processed using the eSASS software system developed by the German eROSITA consortium and proprietary data reduction and analysis software developed by the Russian eROSITA Consortium.

This work includes data products from the NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. The Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged.

This work performed in part under DOE Contract DE-AC02-76SF00515. MG, PM, and RS acknowledge the partial support of this research by grant 21-12-00343 from the Russian Science Foundation. KH has been supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2016R1A5A1013277 and 2020R1A2C1007219), and also financially supported during the research year of Chungbuk National University in 2021.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research has made use of the CIRADA cutout service at URL cutouts.cirada.ca, operated by the Canadian Initiative for Radio Astronomy Data Analysis (CIRADA). CIRADA is funded by a grant from the Canada Foundation for Innovation 2017 Innovation Fund (Project 35999), as well as by the Provinces of Ontario, British Columbia, Alberta, Manitoba and Quebec, in collaboration with the National Research Council of Canada, the US National Radio Astronomy Observatory and Australia’s Commonwealth Scientific and Industrial Research Organisation.

AF received funding from the German Science Foundation DFG, within the Collaborative Research Center SFB1491 ‘Cosmic Interacting Matters – From Source to Signal’. YY thanks the Heising-Simons Foundation for financial support. SR was supported by the Helmholtz Weizmann Research School on Multimessenger Astronomy, funded through the Initiative and Networking Fund of the Helmholtz Association, DESY, the Weizmann Institute, the Humboldt University of Berlin, and the University of Potsdam. ECK acknowledges support from the GREAT research environment funded by Vetenskapsrådet, the Swedish Research Council, under project number 2016-06012, and support from The Wenner-Gren Foundations. MMK acknowledges generous support from the David and Lucille Packard Foundation. This work was supported by the GROWTH project funded by the National Science Foundation under grant no. 1545949.

DATA AVAILABILITY

The data and software to reproduce the main results can be obtained via Zenodo (van Velzen & Stein 2022).

Footnotes

1

While this paper was under review, evidence for a correlation of high-energy neutrinos and a different sample of infrared-selected TDE candidates (Jiang et al. 2021) was presented by Jiang et al. (2023).