© The Authors 2025

1. Introduction

Active galactic nuclei (AGNs) are commonly considered promising candidates for the emission of the high-energy astrophysical neutrinos that have routinely been detected by the IceCube experiment since 2013 (IceCube Collaboration 2013, for a review). Indeed, AGNs can be powerful accelerators of protons, either through relativistic jets or through the accretion disk, hot corona. These particles can then interact with the intense radiation fields present around the active nucleus, thereby producing cascades of particles, including energetic neutrinos. The emission from relativistic jets pointing towards Earth in particular can be very strong due to relativistic boosting (beaming). This is the reason why these ‘oriented’ sources, called blazars, have been extensively investigated in the last ∼10 years, to assess their possible relevance for neutrino astrophysics (Dermer et al. 2014a; Tavecchio et al. 2014; Tavecchio & Ghisellini 2015; Giommi et al. 2020a,b; Giommi & Padovani 2021; Plavin et al. 2020, 2023; Buson et al. 2019, 2022, 2023; Bellenghi et al. 2023).

Blazars can be observed in a wide variety of ‘flavours’ depending on the shape of the spectral energy distribution and/or the optical properties. In particular, blazars in which the synchrotron emission peaks at low (< 10¹⁴ Hz), intermediate (10¹⁴ − 10¹⁵ Hz), or high (> 10¹⁵ Hz) frequencies are called low-, intermediate-, and high-frequency peak sources, respectively (Padovani & Giommi 1995; Abdo et al. 2010). Based on their optical properties, blazars with featureless spectra are instead dubbed BL Lacertae (BL Lac) objects; they are called flat spectrum radio quasars (FSRQs) when strong and broad emission lines, similar to those usually found in non-jetted AGNs, are observed. The distinction between BL Lac and FSRQs could be relevant for neutrino physics since the very existence of a broad-line region (BLR) and a ‘standard’ accretion disk is directly connected to the presence of the radiation field necessary for an efficient neutrino emission (e.g. Dermer et al. 2014b; Murase et al. 2014). For this reason, FSRQs were originally considered the most promising neutrino emitters, although BL Lacs could also be viable candidates if velocity structures are present in the jets (e.g. Tavecchio et al. 2014; Tavecchio & Ghisellini 2015).

Interestingly, the first firm association (> 3σ) between a neutrino event and an extragalactic source, based on the detection of flares temporarily coincident with a neutrino’s arrival and then confirmed by an excess of events in the archival IceCube data, was the source TXS 0506+056 (IceCube Collaboration 2018a,b; Aartsen et al. 2020), an intermediate-frequency peak source that was originally classified as a BL Lac object. However, a detailed analysis of its optical spectrum revealed the presence of a ‘standard’ accretion disk emission (as in FSRQs), whose light was probably swamped by the strong non-thermal continuum coming from the oriented relativistic jet (Padovani et al. 2019). This observation sets TXS 0506+056 apart from ‘true’ BL Lac objects that genuinely lack of an efficient accretion disk and a BLR, probably due to the low accretion rates of their central super massive black holes (SMBHs; e.g. see Hardcastle & Croston 2020), and suggests an emission mechanism that is similar to that considered for FSRQs (Keivani et al. 2018).

The search for positional coincidences between neutrino events (or excesses of events) and extragalactic sources is hampered by the large positional uncertainties of the events. These uncertainties make the number of chance coincidences very high, even for relatively rare objects like blazars.

Hints of positional correlations between neutrino events and BL Lac objects have been found by Padovani et al. (2016) and Giommi et al. (2020b). Other works have suggested a possible association between IceCube events and the most radio-luminous blazars (> 150 mJy), including both FSRQs and BL Lacs (Plavin et al. 2020, 2021, 2023; Buson et al. 2022, 2023).

Even if relativistic jets are natural candidates for neutrino emission, the recent association of a radio-quiet (RQ) AGN (NGC1068) with neutrino events (IceCube Collaboration 2022) suggests that other mechanisms are at work. For instance, Abbasi et al. (2022a) have shown that the acceleration of cosmic rays responsible for the neutrino emission can occur in the accretion disk region of AGNs. The multi-wavelength analysis of NGC 1068 conducted within a p-γ scenario indicates that the most plausible location for the neutrino source is in close proximity to the SMBH, which is embedded in an extremely photon-dense environment (Murase 2022; Padovani et al. 2024). This finding suggests the AGN corona as the likely origin, while star-forming regions, the jet, and the outflow are disfavoured. Yang et al. (2025) find that even in the case of TXS 0506+056, the accretion flow might explain both the long-term steady neutrino emission and the outburst, for which a super Eddington accretion is required. For this reason, the simple positional match of a sample of blazars with neutrino events, even when statistically significant, does not necessarily imply that the emission is due to the relativistic jets, since some of these sources (the FSRQs) also have a powerful accretion disk similar to that observed in RQ quasars.

In this paper we test the hypothesis that FSRQs are responsible, at least in part, for the astrophysical neutrino events detected by IceCube. At the same time, we aim to establish whether this is really due to the presence of an oriented jet or, instead, can be associated with the RQ component of a FSRQ, i.e. the accretion disk and corona. To this end, we started from an optically selected sample of AGNs (that includes both radio-loud and RQ objects), and from this we extracted both the FSRQ and the RQ quasar samples. From an analysis of their correlation with neutrino events, we assessed not only the association of neutrino production with quasars, but also the possible origin site of the neutrinos (jets or disk-corona accretion).

In Sect. 2 we present the quasar data sample and the list of IceCube events that we considered in the analysis. In Sect. 3 we describe the positional quasar–neutrino cross-match and present the results for both FSRQs and RQ quasars. In Sect. 4 we analyse the distribution of the minimum distance between neutrinos and quasars. In Sect. 5 the most likely neutrino emitter candidates are presented individually, and in Sect. 6, we discuss the purity of the sample together with the estimated fraction of high-energy neutrino events produced by FSRQs. Finally, in Sect. 7 we put the results in the context of recent findings presented in the literature.

Throughout the paper we assume a flat Λ cold dark matter cosmology with H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7, and Ω_M = 0.3. Spectral indices are given assuming S_ν ∝ ν^−α.

2. Quasar and neutrino catalogues

2.1. The quasar samples

We started with the sample of type 1 AGNs presented in Shen et al. (2011, hereafter the S11 sample), which is based on the data from the Sloan Digital Sky Survey (SDSS). The sample includes 105 783 spectroscopically confirmed type 1 AGNs, mostly in the northern hemisphere, and contains several pieces of information on the central engine, such as the mass of the SMBH and its accretion rate, which can be used to further investigate potential neutrino emitters. For our analysis we only considered the sky area within 120 ° < RA < 260° and 0° Dec 60°, where the coverage is uniform (Fig. 1); this area includes a total of 89 437 quasars.

Fig. 1. Sky distribution of the data. The small grey dots represent the positions of 105 783 quasars from the S11 catalogue. The small blue crosses indicate the 1218 CLASH FSRQs. Black circles mark the 33 IceCAT-1 events in the northern sky with the required positional accuracy. Circles with a red centre represent the 15 neutrinos located within our area of interest uniformly covered by the CLASH objects. A red square around a circle indicates that a CLASH FSRQ is included in the nominal error region.

For the radio information we used the Cosmic Lens All Sky Survey (CLASS; Myers et al. 2003; Browne et al. 2003), which homogeneously covers all the extragalactic northern sky. The CLASS catalogue contains ∼11 000 flat spectrum radio sources selected at 5 GHz down to a well-defined flux density limit (30 mJy). Since blazars are flat spectrum radio sources, this survey represents the best starting point for our selection. The survey is also deep enough to select most of the blazars with an optical magnitude brighter than 20 (which roughly corresponds to the limit of the S11 sample). Thanks to the excellent positional accuracy of the survey (< 1″) the optical/radio association can be performed with high reliability (Fig. 1). The cross-correlation between the S11 sample and CLASS produces 1218 sources (1122 in the region of interest), using 1″ positional tolerance. We call this sample CLASH, and throughout the rest of the paper, unless otherwise specified, flux densities are measured at 5 GHz.

For comparison, we also built a sample of RQ quasars, defined as those sources from S11 not detected in CLASS and located more than 5″ from any source in the FIRST catalogue. FIRST is a radio survey at 1.4 GHz with flux density limit of ∼1 mJy, significantly deeper than CLASS and with a good positional accuracy (∼2″). Given the typical magnitudes of the Shen et al. (2011) quasars, the non-detection of a source in the FIRST survey implies a radio-loudness¹ below 10 i.e. in the typical RQ regime (Kellermann et al. 1989). Therefore, the RQ quasar sample should mostly contain non-jetted objects. This control sample comprises 73 339 RQ objects.

As shown in Fig. 2, the two populations are indistinguishable in terms of redshift, mass, bolometric luminosity and Eddington ratio. The only difference is the presence of a strong radio emission, which we assume is produced by a relativistic jet oriented towards our line of sight.

Fig. 2. Comparison of the properties of CLASH and RQ quasars. The distributions of the redshift, mass, bolometric luminosity, and Eddington ratio of the RQ sample (black) and CLASH blazars (blue) are shown. The values of the possible neutrino emitters are given in red (Table 3).

2.2. The IceCube neutrino events

For our study we used the IceCube Event Catalog of Alert Tracks (IceCat-1; Abbasi et al. 2023b), which includes 267 non-vetoed track-like neutrino events likely produced by astrophysical sources between 2011 and 2021. Since the positional uncertainty is extremely large in a significant fraction of the events, we considered only those events with the best positional accuracy:

$$\begin{array}{cc}\hfill ({\mathrm{ERR}}_{\mathrm{RA}}^{-}+{\mathrm{ERR}}_{\mathrm{RA}}^{+})\cos (\mathrm{Dec})& <2.0^{°}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mathrm{ERR}}_{\mathrm{Dec}}^{-}+{\mathrm{ERR}}_{\mathrm{Dec}}^{+}& <2.0^{°}.\hfill \end{array}$$(2)

We stress that these requirements were frozen at the beginning of the analysis, and no optimization has been performed. In the region of interest defined above (Fig. 1), these constraints produce a list of 15 events (Table 1), with IC121103A, IC190503, and IC200117A excluded because too close to the border. Following the approach of Plavin et al. (2023), also discussed in Abbasi et al. (2023a), the 90% error region is assumed to be encircled by the 4 elliptical segments with the RA and Dec. errors as semi-axes. The energy distribution of the 15 events is shown in Fig. 3: the selection on the positional accuracy includes events belonging to the bulk of the population detected in the northern sky with energies in the 100–500 TeV interval (Table 1).

Fig. 3. Energy distribution of the 15 selected events (red) compared to the energy distributions of all 267 Icecat-1 events (grey) and those observed in the northern sky (hatched grey).

3. Quasar–neutrino associations

3.1. CLASH FSRQ–neutrino association

Considering the list of 15 IceCat-1 events and the 1122 CLASH FSRQ we find that 9 CLASH sources fall within the error boxes of 7 IceCube events (Table 1). Varying the flux density limit of the CLASH sample, the number of objects varies from 1122 with a limit of 30 mJy to 270, when the flux density limit is 200 mJy. The number of neutrino events, with at least one CLASH quasar within the 90% error region varies from 7 to 3, while the number of CLASH included in the neutrino 90% confidence error region varies from 9 to 3 (Table 2).

Multiple matches clearly indicate that the probability of chance coincidence is not negligible. We estimated the expected number of chance coincidences by randomly shifting quasar positions within a circular annulus of radius 3–10° around their original locations. This procedure was repeated 10⁵ times for different radio flux density limits of the FSRQ sample. The statistical significance of the observed correlation (p-value) was calculated by counting the realizations in which the number of neutrinos including at least one CLASH FSRQ exceeded that observed in the actual dataset. The output of this Monte Carlo simulation at different radio flux density limits are reported in Fig. 4 and Table 2.

Fig. 4. p-values for the neutrino–CLASH associations (red) and neutrino–S11-RQ quasar associations (grey). The solid red line, plotted for different flux density limits, represents the probability of randomly finding an equal or greater number of neutrino events with at least one CLASH FSRQ within the error box. The dotted and dashed lines indicate the probabilities that the distribution of minimum distances between neutrinos and CLASH sources arises from a random distribution (KS test) for all 15 events and for the high-declination sample, respectively. The p-values for 100 RQ samples, of the same size as the FSRQ sample at the corresponding radio flux limit, are displayed in grey. The median value, along with the 90th and 99th percentiles, is calculated from 10 000 samples, which are not shown for clarity. We plot the p-values for the most luminous and most massive quasars of the S11 sample in green and blue, respectively.

These values show a minimum located at a flux density limit of ∼100 mJy (Fig. 4). At this flux density limit, the number of CLASH FSRQs is 497, with 7 of them in the error box of 6 neutrino events. The average number of associations of an equal number of random positions with the 15 neutrino error box is ∼1.3. We find that in 52 cases out of 1e5 an equal or larger number neutrino have a CLASH in the error box. Considering the number of different trials involved in this analysis (6 with different flux density limits), we estimate that the final probability corresponds to a significance of ∼2.7σ. We obtain identical results when repeating the simulation by randomly shifting the positions of the neutrinos instead of the quasars.

3.2. RQ quasar–neutrino association

Among the 77339 RQ quasars, 178 fall within the 15 neutrino error regions. Since these numbers do not allow the assessment of any possible significant correlation, we followed two different approaches. First, we compared the RQ quasars to the FSRQ sample with similar statistic, using the 77 339 RQ sources to construct 10 000 different random samples of 1122 RQ objects each. For each sample, we repeated the same analysis performed on the CLASH FSRQs, extracting subsamples of the same size as the CLASH sample at different flux density limits, ranging from 1122 to 270 objects. For each of these samples, we counted the matches with the 15 neutrino positions and estimated the expected number of chance coincidences again by randomly shifting the positions. We then repeated the same analysis on the 1,122 most massive and the 1122 most luminous quasars, regardless of their radio flux.

The p-value results are shown in Fig. 4. We find average p-values of ∼0.7 with the 99% of the samples having p-values higher than 0.02, meaning that the number of RQ quasars present in the neutrino error box is fully consistent with what expected from chance coincidence.

The absence of a correlation between RQ quasars and neutrino events, combined with the correlation found for the FSRQ subsample, strongly suggests that, if FSRQs are confirmed as neutrino sources, the production mechanism must be tied to the physics of the jet rather than the accretion flow.

4. Neutrinos and the closest CLASH FSRQs

In the previous section we assessed the correlation between neutrinos and CLASH objects by counting the number of neutrino events that include FSRQs within their 90% confidence error regions. In this section we present a different approach, analysing the angular distances between each of the 15 considered neutrino events and the nearest CLASH source. Figure 5 presents the minimum distances between the 15 neutrino positions and CLASH sources as a function of declination, considering only objects with a flux density above 100 mJy (497 sources). The figure also shows, for comparison, the mean and standard deviation of the distribution of minimum distances obtained from random control samples. Specifically, we generated 10 000 control samples of 497 elements, randomly extracted from the RQ population, which uniformly covers the same area as the CLASH sample. As an additional check, we produced 10 000 further samples by randomizing the positions within a circular annulus of radius 3–10° around the original locations of the CLASH and RQ quasars, finding perfectly consistent numbers.

Fig. 5. Distances between the 15 neutrino events and the closest CLASH FSRQs as a function of the declination (red circles). A square around the circle indicates that the CLASH FSRQ is included in the nominal error region. For each event, in grey, we plot the mean and standard deviation of the distribution of the minimum distances of a quasar sample of the same size with randomized positions.

A clear trend emerges when comparing the subsamples at low and high declinations. At low declinations, the distances between seven out of the eight neutrino events and the nearest CLASH source are consistent with random expectations. In contrast, at higher declinations, all seven neutrino events have a nearby CLASH bright object within ≲0.7°, significantly closer than expected from random distributions.

We applied a Kolmogorov-Smirnov (KS) test to evaluate whether the CLASH-neutrino distance distribution differs significantly from a random distribution. The results, for both the full sample of 15 neutrino events and the high-declination subsample (7 events), are presented in Fig. 4 and summarized in Table 2 for different radio flux density thresholds. While the difference is moderately significant for the full sample, the high-declination subsample shows a much stronger deviation from randomness. The minimum uncorrected p-value at flux limit of 100 mJy is 4 × 10⁻⁶, corresponding to a significance of 4.0σ after accounting for the six trials (Fig. 4).

The dependence of the distance between the closest CLASH-FSRQ and neutrino events on declination may have a technical explanation. The performance of the IceCube telescope, particularly in terms of background levels and effective area, varies significantly with declination (Abbasi et al. 2023b), which corresponds to the (opposite of the) neutrino incident direction. Although we are not in a position to fully assess the event detection error budget, our results suggest that the positional accuracy and purity level of IceCat-1 catalogue events are higher at declinations above 20°.

As the number of events with good positional accuracy in IceCat-1 within our region of interest is limited, in order to check this result, we expanded our sample. We considered the entire CLASS catalogue, regardless of identification with an S11 quasar. To be homogeneous with the S11 catalogue, we considered only sources with magnitudes of < 20 in at least one of the SDSS band. With this limit at flux densities of > 100 mJy, the CLASS catalogue includes 1790 objects, of which 1052 are classified as FSRQs in the literature; the remaining objects are unclassified.

This extension allowed us to use all the 33 events in the northern hemisphere with good positional accuracy (as defined earlier). Additionally, we included the 18 neutrino events that meet the same positional accuracy criteria and occurred from January 2022 to December 2024, as reported in the GCN catalogue². The distances between these 51 events and the closest CLASS source are reported in Fig. 7. Among the four objects added to the high-declination sample, three have a minimum distance ≲0.4°, while one is at 1.8°. We note that one of the three at a short distance is GB6J0740+28, a CLASH object in the error region of IC200117A, which was excluded from the sample of 15 events because it is too close to the border of our region of interest (Sect. 2.2). Therefore, the pattern of high-declination events being closer to a bright radio source seems to also be present in this expanded catalogue.

5. High-declination sample

Focusing on the original sample of 15 events, we find that, among the seven neutrinos detected at declinations above 20°, five have a bright CLASH FSRQ (flux >  100 mJy) within their 90% confidence error region, while all seven have at least one CLASH FSRQ within 0.7° (Fig. 5). We note that the two high-declination events without a bright CLASH source within their nominal 90% confidence error region would include such a source if a systematic error of 0.1° were (linearly) added to their nominal positional uncertainties (Fig. 6).

Fig. 6. Maps of the seven CLASH objects identified as possible neutrino sources. Black circles show nominal neutrino positions. Dashed lines show the 90% confidence error regions, and the solid lines show the elliptical approximation. The CLASH FSRQs are indicated by red crosses. Green squares show BZCat blazars, and purple ellipses the Fermi catalogue element 90% confidence error regions. The small blue crosses are the S11 quasar positions.

Fig. 7. Same as Fig. 5 but for a larger sample of neutrino events. The radio counterparts in this case are taken from the CLASS catalogue, regardless of identification with an S11 quasar.

In this section we highlight the main characteristics of these seven CLASH sources. All the information presented here is derived from public databases and catalogues. The seven FSRQs have redshifts in the 0.5–2.3 range and optical AB magnitudes between 17 and 19. X-ray fluxes in the 0.5–2.0 keV band vary between 1 and 5 × 10⁻¹³ erg s⁻¹ cm⁻², which are typical values for FSRQs with radio fluxes > 100 mJy. The SMBH masses range between 10^8.8 and 10^9.7 solar masses and are accreting at rates estimated between 0.1 and 0.8 in terms of Eddington ratios (Table 3). We investigated whether these numbers are consistent with the characteristics of the S11 RQ quasar population. The analysis is hampered by the small sample size and no significant difference was detected, with K-S test probabilities exceeding 10% in all three cases.

From a visual inspection of the spectral energy distribution of these seven objects, determined using the FIRMAMENTO³ tool (Giommi 2025), we find that they all have a usual blazar shape, with a synchrotron peak falling at 10¹³–10¹⁴ Hz, i.e. in the typical range observed in FSRQs. Only one object listed in Table 3 was detected in gamma rays by Fermi; it is listed as positionally associated with a source in the 4FGL-DR4 catalogue (Abdollahi et al. 2020; Ballet et al. 2023). A second object (GB6J1125+2610) is indicated only as a possible, low-confidence, association. Therefore, only a small fraction of the FSRQs that are likely associated with a neutrino event seem to be clearly detected in gamma rays. Neutrinos produced via photo-pion interactions are expected to be accompanied by an equivalent emission of γ-ray photons. However, the neutrino production site is likely optically thick to high-energy photons (e.g. Böttcher 2019), and therefore, the gamma rays produced co-spatially with neutrinos are expected to be strongly absorbed and not detectable. In this scenario, the observed γ-ray emission can be attributed to the leptonic component through the inverse Compton mechanism. FSRQs with a low frequency peak (< 10¹³ − 10¹⁴ Hz) usually have the high-energy hump peaking at energies ≤100 MeV, i.e. below the sensitivity of Fermi. These blazars could appear weak in Fermi, even if their luminosity at the peak is high. Notably, the detection of one (possibly two) sources is consistent with the detection rate of bright FSRQs (flux > 100 mJy) by Fermi, which stands at approximately 40%.

One of these objects (GB6J1125+2610, also known as TXS 1123+264) has been previously indicated as a possible neutrino emitter by Plavin et al. (2020). Interestingly, a study of 20 years of archival very long-baseline interferometry (VLBI) data detected a significant brightening of the core after the neutrino event (Kőmíves et al. 2024).

Moreover, we compared the positions of the seven FSRQ candidates with the 9-year IceCube sky map (IceCube Collaboration 2022), which provides the probability that the observed neutrino flux is of non-astrophysical origin, mapped across the entire sky on a 0.2 ° ×0.2° grid. Interestingly, we find that four out of the seven FSRQ candidates lie within 1° of a significant peak (local p-value < 4 × 10⁻⁴) in the sky map To calculate the number of associations expected by chance, we exploited the uniform distribution of RQ quasar positions over the region by extracting 10 000 samples of seven positions each. We find a mean value of 0.12. The most intriguing case is GB6J1350+233 (CRATESJ1350+23), located approximately 0.8° from the fifth most significant position on the IceCube scan (IceCube Collaboration 2022, Table S2). Although the sky map and the IceCat-1 catalogue are not independent datasets, this correlation might suggest that, with a high probability, some of these seven FSRQs produce detectable and measurable neutrino fluxes.

6. Bronze and gold events and the fraction of neutrinos produced by FSRQs

In the IceCat-1 catalogue, neutrino events are classified as either bronze or gold, based on their estimated probability of being genuinely astrophysical in origin. Bronze events are assigned a 30% probability, while gold events have a 50% probability (Abbasi et al. 2023b). In the sample of 15 neutrino events considered in our analysis, six are flagged as gold and nine as bronze, leading to an expected total of 5.7 real astrophysical neutrino events. If we define an association as a positional match within 0.7°, 8 out of the 15 neutrinos have at least one associated CLASH FSRQ, with one event (IC180125A) having two associations. Using our procedure to generate control samples with randomized quasar positions, we estimate that only 1.5 such associations are expected by chance. This implies 7–8 genuine associations, which is consistent with the expected 5.7 astrophysical, at least under the hypothesis that a significant fraction of them are produced by FSRQs. However, focusing only on the high-declination subsample yields a different picture. As noted, all seven neutrino events located at declinations above 20° are associated with a CLASH FSRQ, compared to only 0.8 associations expected by chance. Among these, five events are bronze and two are gold, implying an expected 2.5 real astrophysical events, inconsistent with the observed 7 associations (probability ≲ 1%). As discussed earlier, the observed difference between the high- and low-declination samples may be due to unaccounted systematic uncertainties, affecting either the positional accuracy or the purity of the IceCat-1 catalogue. The fact that all seven high-declination events are associated with a nearby FSRQ suggests that this subsample has a higher astrophysical purity.

Assuming that all the 7 neutrino are real, this would imply that the fraction of neutrino produced by FSRQs is > 60% at 2σ confidence. Such a high fraction of high-energy neutrino events originating from FSRQs appears to be in tension with the much lower values reported in other studies (< 20%, Palladino et al. 2020; Oikonomou 2022, and references therein). The first argument in those studies is based on the expected electromagnetic counterpart, primarily in the gamma-ray band. These works typically estimated upper limits on the total neutrino flux from blazars by assuming a fixed neutrino/γ-ray flux ratio. However, as discussed in the previous section, many FSRQs can appear faint or even undetectable by Fermi, as their high-energy spectral peak can fall below the Fermi LAT sensitivity range. As a result, the contribution of such sources to the neutrino flux could be significantly underestimated in gamma-ray-based analyses.

A second argument comes from the lack of significant neutrino point sources, which, combined with the observation of a diffuse astrophysical neutrino flux, implies that the dominant sources should form a dense population of weak emitters, rather than a few luminous ones. In a recent analysis, Abbasi et al. (2023c) found that, assuming steady emission and evolution following the star formation rate, a minimum local source density of 7 × 10⁻⁹ Mpc⁻³ is required to match the observed neutrino flux. FSRQs, with a local density of only 5.6 × 10⁻¹¹ Mpc⁻³, as inferred from the Fermi catalogue, are too rare to contribute significantly. However, this conclusion is highly dependent on the assumed source evolution. In the same study (Abbasi et al. 2023c), the authors note that adopting the evolution of the AGN luminosity function would reduce the required source density by nearly two orders of magnitude, down to 6 × 10⁻¹¹ Mpc⁻³, much closer to the FSRQ density.

In our analysis we find that the neutrino-associated FSRQs have radio luminosities at 1.4 GHz in the range of 10⁴²–10⁴⁴ erg s⁻¹. The space density and redshift evolution of this population have been well characterized in previous works (Mao et al. 2017; Caccianiga et al. 2019; Diana et al. 2022). Their local space density is estimated to be around 4 × 10⁻¹⁰ Mpc⁻³, with a strong positive evolution peaking at z ∼ 1.9, where the density reaches 2.5 × 10⁻⁸ Mpc⁻³. This evolution is significantly faster than that of the star formation rate and more closely resembles that of AGNs. Taking this into account, our findings are consistent with both the absence of neutrino multiplets and the observed diffuse neutrino flux and thus reconcile the high apparent association rate of FSRQs with the broader neutrino source population. A similar result has been found by Neronov & Semikoz (2020).

Abbasi et al. (2022b), using a stacking analysis of 137 FSRQs selected in the 30–100 MeV band, found that their contribution to the unresolved background is ≲1%. At first glance, this result may seem to be in tension with our findings, since some level of soft γ-ray/hard X-ray emission is inevitably expected from FSRQs if they are indeed neutrino emitters (see the discussion in Abbasi et al. (2022b). However, the expected contribution of these very bright and rare objects to the unresolved neutrino flux is not significant. Furthermore, the hadronic contribution to the soft γ-ray flux is difficult to quantify, while the inverse Compton emission is certainly contributing in this band. It is therefore plausible that the 137 FSRQs bright in the soft γ-ray band considered by Abbasi et al. (2022b) are simply those with the strongest Compton dominance, and thus not necessarily the most promising neutrino sources. We do expect that the FSRQs discussed in this work emit in the soft γ-ray band, but it is likely that the sensitivity of Fermi in this energy range is insufficient to detect them.

7. Discussion

Several studies aimed at the search for positional correlations between IceCube events and extragalactic sources were carried out in the last few years. They all adopted different approaches and started from different neutrino events and AGN catalogues. In the following sections we summarize the main results coming from the studies.

7.1. Neutrino events

Concerning neutrino events, many authors (including us) considered the list of track-like events with relatively good positional accuracies (Giommi et al. 2020a; Plavin et al. 2020, 2023). Early searches used cascade-type events (Padovani et al. 2016), which are characterized by larger positional uncertainties but have a more likely astrophysical origin. Apart from the list of single neutrino events, the IceCube collaboration has also provided an all-sky map of possible neutrino excesses (Aartsen et al. 2017) in the form of a p-value, for each sky position, that gives the level of clustering in the neutrino data, i.e. a measure of the significance of neutrino events being uniformly distributed. The map is based on more than 700 000 events taken during 7 years of data collection. This map was used by the IceCube collaboration (Aartsen et al. 2017) and by other authors (Plavin et al. 2021; Buson et al. 2023; Bellenghi et al. 2023) to search for astrophysical counterparts. More recent versions of the map, based on 10 years of data, were produced and used by Aartsen et al. (2020) and Bellenghi et al. (2023).

7.2. FSRQs versus BL Lacs

As far as the blazar catalogues are concerned, different possibilities were considered. A common choice was the selection of gamma-ray blazars, detected in one of the available Fermi catalogues (Aartsen et al. 2017, 2020) or blazars (typically BL Lacs) that are promising candidates to be detected at very high-energy gamma rays (Giommi et al. 2020b). This choice was motivated by the fact that, in the hadronic model, a neutrino emission should be accompanied by the emission of a gamma-ray photon. However, not only the gamma-ray photon can be absorbed, right after its emission or during its travel towards Earth (by the diffuse extragalactic background light), but a gamma-ray emission can be also produced by other mechanisms (e.g. the inverse Compton scattering) without the production of a neutrino (leptonic model). It is possible that both leptonic and hadronic mechanisms are at work in the jets of blazars, and therefore, a selection based on gamma-ray strength may not be the optimal choice. A gamma-ray ‘agnostic’ approach was adopted by Buson et al. (2022, 2023) and Plavin et al. (2020, 2021, 2023) who considered, respectively, all the blazars included in the Roma-BZCat catalogue (5BZCat; Massaro et al. 2015) and the radio brightest blazars, based on the Very Long Baseline Array calibrator surveys. In all these works, blazars are analysed as a unique class, without distinguishing between BL Lac objects and FSRQs.

Overall, the most significant (≥3σ post-trial) correlations found so far are those with the BZCat blazars (Buson et al. 2022) and with radio bright (> 150 mJy, Plavin et al. 2020, 2021) blazars. In all these works, the FSRQs that seem to correlate with the neutrino events are radio bright, above ∼200 mJy. In the ‘dissection’ analysis discussed in Giommi et al. (2020a), a > 3σ (post-trials) correlation is found only with blazars (including both BL Lacs and FSRQs) with the synchrotron component peaking at high frequencies (> 10¹⁴ Hz) while no significant correlation is found with low-frequency peaked blazars. Fermi-bright blazars do not seem to provide significant correlations except when using the most recent (10 years) map (Aartsen et al. 2020). The correlation with radio bright blazars from the RFC catalogue (Petrov & Kovalev 2025) was confirmed at 3.3σ also in the southern hemisphere by Bellenghi et al. (2023) using the 7-year map although nothing significant was found using the revised map based on 10 years of data. This could be in part due to the lower sensitivity achieved by their analysis compared to Aartsen et al. (2020) caused by the lack of an accurate detector response matrices publicly available (Bellenghi et al. 2023).

In the work presented here, we obtain a significant correlation (∼4σ) based on a sample of radio bright (> 100 mJy) FSRQs with declinations above 20°. The type of selection is similar to that discussed in Plavin et al. (2020) although we focused only on the FSRQs, excluding BL Lac objects. Moreover, we restricted the analysis to the sky area covered by SDSS in order to have an almost complete and uniform spectral identification of the radio bright sources. It is interesting to note that the four sources in the Plavin et al. (2020) that are indicated as most likely counterparts of the neutrino events are all classified as FSRQs, even if the starting sample includes both FSRQs and BL Lac objects. These four FSRQs are not included in our sample since they fall in an area of sky that we are not covering. Our result is corroborated by the observation that at high declination (Dec > 20°) the significance of the correlation is higher. Based on these results, it seems that the radio-brightest FSRQs, mostly of low-frequency peak type, give a high signal in the correlation analyses carried so far.

We note that the ‘dissection’ analysis performed by Giommi et al. (2020a) did not find any significant correlation when considering low-frequency peaked blazars. This could be due to the fact that the objects considered in that work include both BL Lacs and FSRQs and, more importantly, radio faint sources. For a direct comparison with our findings, this kind of analysis should only be carried out on bright FSRQs (> 100 mJy).

We finally note that TXS0506+056 is formally classified as a BL Lac object and, therefore, a priori excluded in the analysis presented here. The fact that broad emission lines, although present (Padovani et al. 2019), were not visible in the discovery spectra is related to the fact that the synchrotron component peaks at a relatively high frequency peak (10¹⁴–10¹⁵ Hz, Padovani et al. 2018), i.e. in the visible range, which is higher compared to ‘normal’ FSRQs, where this component typically peaks in the far infra-red band. Because of this peculiar characteristic, the non-thermal emission from the jet is able to swamp the accretion disk and BLR components, thus hiding the emission lines. These objects, also called ‘masquerading BL Lacs’ (Giommi et al. 2013), represent a rare tail of the entire FSRQ population. Distinguishing these sources from ‘true’ BL Lacs is not simple, but their inclusion in a future analysis might increase the significance of our results.

7.3. Jet versus accretion

As mentioned in the Introduction, neutrino production mechanisms linked to matter accretion onto the SMBH have been proposed in the literature for both NGC 1068 and TXS 0506+056. However, the results presented in this paper suggest that the neutrino production mechanism in these FSRQs is not associated with accretion. The correlation between FSRQ positions and neutrino events, in contrast with the lack of correlation of samples of quasars with similar optical properties but no radio emission, strongly indicates that neutrino emission originates in relativistic jets. To our knowledge, it is the first time that such a direct comparison between radio-loud and RQ quasars, in the search for counterparts of neutrino events, is carried out.

We note that our result is based solely on the positional correlation between neutrino events and radio jet emission. Other studies have found evidence of temporal correlations between radio activity and neutrino detections (Plavin et al. 2020). Since the radio flux densities used in our analysis were measured at epochs completely independent of the neutrino events, our approach does not allow us to test for a direct temporal connection between radio emission and neutrino production. In the lepto-hadronic scenario, which appears to provide the best match to observational data, the radio emission is primarily produced by synchrotron radiation from the leptonic component (e.g. Petropoulou et al. 2020) and is not directly associated with the pγ interactions responsible for neutrino generation. The high radio flux levels observed in our analysis can be interpreted as a selection effect, with neutrinos being more likely detected from jets that are either intrinsically more powerful or strongly Doppler boosted (Plavin et al. 2023).

8. Summary and conclusions

Starting from an optically selected sample of AGNs (Shen et al. 2011) and a well-defined, flux-density-limited radio survey at 5 GHz in the northern hemisphere (CLASS), we identified a sample of 1122 FSRQs, referred to as the CLASH sample, from the larger RQ quasar population with identical optical properties. We then performed a simple positional cross-correlation between both samples and a set of 15 IceCube events extracted from IceCat-1, selecting events with relatively well-constrained positions within the SDSS/CLASS survey area. Using a Monte Carlo randomization of positions, we estimated the p-values of the observed correlations. Our analysis of the CLASH sample at six different radio flux density limits reveals that the most significant correlation occurs at a flux density limit of 100 mJy, corresponding to a post-trial significance of 2.7σ. In contrast, applying the same test to the RQ quasar population yielded significantly higher p-values that were fully consistent with a random distribution.

Analysing the distribution of minimum distances between neutrino events and CLASH sources, we observe a clear trend: at declinations above 20°, these distances are systematically smaller than what would be expected by chance, whereas closer to the celestial equator, they tend to be larger. A KS test confirms this trend. While the minimum distance distribution for the full sample differs from a random distribution at a moderate significance level of 2.2σ, restricting the analysis to high-declination events increases the statistical significance to 4.0σ. In contrast, we find that the minimum distance distribution for the RQ quasar sample is indistinguishable from a random sample. This reinforces the observed correlation between neutrinos and FSRQs.

These findings indicate that a majority (> 60%) of the high-energy astrophysical neutrinos detected by IceCube originate from blazars, particularly FSRQs. Their production likely occurs in relativistic jets rather than in the RQ regions, such as the accretion disk or corona. Our expectation is that most neutrino events that will be detected at high declinations by IceCube, particularly those with relatively good positional accuracies, will be located close to a luminous source with a flat radio spectrum and broad emission lines in the optical.

Acknowledgments

The authors thank Paolo Giommi and Lorenzo Caccianiga for valuable discussions and suggestions. We acknowledge the financial support from INAF under the projects “Quasar jets in the early Universe” (Ricerca Fondamentale 2022) and “Testing the obscuration in the early Universe” (Ricerca Fondamentale 2023).

References

All Tables

All Figures

Fig. 1. Sky distribution of the data. The small grey dots represent the positions of 105 783 quasars from the S11 catalogue. The small blue crosses indicate the 1218 CLASH FSRQs. Black circles mark the 33 IceCAT-1 events in the northern sky with the required positional accuracy. Circles with a red centre represent the 15 neutrinos located within our area of interest uniformly covered by the CLASH objects. A red square around a circle indicates that a CLASH FSRQ is included in the nominal error region.

In the text

	Fig. 2. Comparison of the properties of CLASH and RQ quasars. The distributions of the redshift, mass, bolometric luminosity, and Eddington ratio of the RQ sample (black) and CLASH blazars (blue) are shown. The values of the possible neutrino emitters are given in red (Table 3).
In the text

	Fig. 3. Energy distribution of the 15 selected events (red) compared to the energy distributions of all 267 Icecat-1 events (grey) and those observed in the northern sky (hatched grey).
In the text

Fig. 4. p-values for the neutrino–CLASH associations (red) and neutrino–S11-RQ quasar associations (grey). The solid red line, plotted for different flux density limits, represents the probability of randomly finding an equal or greater number of neutrino events with at least one CLASH FSRQ within the error box. The dotted and dashed lines indicate the probabilities that the distribution of minimum distances between neutrinos and CLASH sources arises from a random distribution (KS test) for all 15 events and for the high-declination sample, respectively. The p-values for 100 RQ samples, of the same size as the FSRQ sample at the corresponding radio flux limit, are displayed in grey. The median value, along with the 90th and 99th percentiles, is calculated from 10 000 samples, which are not shown for clarity. We plot the p-values for the most luminous and most massive quasars of the S11 sample in green and blue, respectively.

In the text

	Fig. 5. Distances between the 15 neutrino events and the closest CLASH FSRQs as a function of the declination (red circles). A square around the circle indicates that the CLASH FSRQ is included in the nominal error region. For each event, in grey, we plot the mean and standard deviation of the distribution of the minimum distances of a quasar sample of the same size with randomized positions.
In the text

Fig. 6. Maps of the seven CLASH objects identified as possible neutrino sources. Black circles show nominal neutrino positions. Dashed lines show the 90% confidence error regions, and the solid lines show the elliptical approximation. The CLASH FSRQs are indicated by red crosses. Green squares show BZCat blazars, and purple ellipses the Fermi catalogue element 90% confidence error regions. The small blue crosses are the S11 quasar positions.

In the text

	Fig. 7. Same as Fig. 5 but for a larger sample of neutrino events. The radio counterparts in this case are taken from the CLASS catalogue, regardless of identification with an S11 quasar.
In the text