© I. Liodakis et al. 2022

1. Introduction

Throughout the history of modern astronomy, the introduction of new theories and new technologies has led to giant leaps in our understanding not only of astronomy, but also of physics and cosmology. The successive openings of the radio, X-ray, γ-ray, and most recently the gravitational wave windows, have revolutionized our view of the Universe. The opening of the high-energy ∼PeV neutrino window is likely among the last observational frontiers. It is reasonable to expect that opening this window will lead to fundamental changes in our understanding, as has happened in all previous cases. The start of operations of IceCube in 2010 and the detection of high-energy astrophysical neutrinos (IceCube Collaboration 2013; Aartsen et al. 2014) was a significant step towards this goal. Identifying sources of neutrinos is the next obvious, yet paramount, step forward.

In the one proven case of multi-messenger astrophysics, that is that of gravitational waves, each potential gravitational wave source has a distinctive gravitational wave signature. Therefore we know what we are looking at and what we are looking for. In the case of multi-messenger neutrino astrophysics there are no distinctive neutrino signatures. This is exacerbated by the fact that at high energies we may only receive one neutrino from each source, which is not sufficient to provide a neutrino signature. To make things worse, there is a huge background of atmospheric muons (∼10¹¹ yr⁻¹) and atmospheric neutrinos (∼10⁵ yr⁻¹) against which the high-energy astrophysical neutrinos (∼10  yr⁻¹) have to be detected. Nevertheless, in 2017 IceCube observed a high-energy neutrino event designated as IceCube-170922A. The best-fit reconstructed direction was at 0.1° from the sky position of the blazar [an active galactic nucleus (AGN) with a jet pointed towards the observer] TXS 0506+056, which was flaring in radio, in X-rays and in GeV and TeV γ-rays (Tanaka et al. 2017; Mirzoyan 2017; Kun et al. 2019). The positional and temporal coincidences were significant at the 3σ level, suggesting the detection of the first counterpart to a high-energy extragalactic neutrino (IceCube Collaboration 2018). There is one other radio bright jetted-AGN within 2 degrees of the IceCube position, and the number of sources is expected to scale as S^−3/2.

In relativistic jets, there is now evidence that electromagnetic (EM) emission can originate in the outer surface of the jet, enhanced by interaction with the interstellar medium (ISM, Kim et al. 2018; Janssen et al. 2021). In this situation, one could expect entrainment of ISM material, and hence of hadrons, that would give rise to a hadronic emission channel in relativistic jets. Thus relativistic jets in AGN could produce neutrinos and it is concerning that this association remains in doubt after over a decade of observations by IceCube and the detection of dozens of high-energy astrophysical neutrinos. As confidence in the existence of astrophysical high-energy neutrinos grows through the increasingly sophisticated removal of the contaminating atmospheric muons and neutrinos, if no EM counterparts are found, the possibility of dark matter sources and of physics beyond the standard model gains strength. These are compelling arguments for mounting the most powerful attacks possible for the identification of the EM sources of astrophysical neutrinos. This paper has been written partly in order to draw attention to this strategic deficit in the pursuit of multi-messenger neutrino astrophysics.

Discussions of high-energy (TeV–PeV) neutrinos in relativistic jets have involved pion production through proton-proton (e.g., Loeb & Waxman 2006) and proton-photon interactions (e.g., Stecker et al. 1991; Mannheim 1995; Halzen & Zas 1997). This requires proton acceleration to energies typically an order of magnitude larger than the observed neutrino energies. Interestingly, these processes predict an initial flavor ratio of ν_e:ν_μ:ν_τ of 1:2:0 for pp collisions that produce roughly equal number of π+ and π−. For pγ interactions close to the Δ resonance, the relative ratio can be 1:1:0 (Anchordoqui et al. 2005). The neutrinos then mix, in an energy-dependent manner, on the journey to Earth. The mixing ratios could be confidently studied by an expanded IceCube facility such as IceCube-Gen2 (The IceCube-Gen2 Collaboration 2021).

Efficient particle acceleration in blazar jets is quite natural. The most luminous flat spectrum radio quasars (FSRQ) near their source can generate up to ∼1 ZV = 10²¹ V of electromotive force and, for this reason, quasar jets have commonly been suspected as the source of ultra high energy cosmic rays with energy up to ∼10²⁰ eV (e.g., Globus et al. 2017; Mbarek & Caprioli 2021). It is therefore quite reasonable that protons, given the longer cooling timescales than leptons, can be accelerated to high enough energies to account for the observed astrophysical neutrinos. However, the charged pions produced in these interactions also create pairs which will radiate γ-rays either through inverse-Compton scattering or synchrotron radiation. Neutral pions create γ-rays directly. Following the detection of very-high-energy γ-rays and the IceCube-170922A event from the direction of the jetted-AGN TXS 0506+056, various theoretical models have been applied to explain the observed association of neutrinos and blazar flares (see e.g., Ansoldi et al. 2018; Sahakyan 2018; Murase et al. 2018; Keivani et al. 2018; Zhang et al. 2019; Petropoulou et al. 2020) as well as neutrinos and steady emission processes of blazar jets (see e.g., Murase et al. 2018; Rodrigues et al. 2019). Many of these papers show that it is not easy to explain the multi-messenger data of TXS 0506+056. First, as suggested for example by Keivani et al. (2018) and Zhang et al. (2019), a physically consistent picture can only be achieved with a hybrid leptonic scenario, with γ-rays produced by external inverse-Compton processes and high-energy neutrinos produced by photopion production with the external radiation field. The best-fit models, however, typically produce neutrinos peaked at higher energy but with lower flux compared to observations. Alternatively, a hadronuclear interaction channel has also been investigated (see e.g., Liu et al. 2019). However, producing the observed neutrino flux through proton-proton interactions requires an extreme jet power for accelerating relativistic protons that significantly exceeds the Eddington luminosity (see also, Liodakis & Petropoulou 2020). In addition single-zone models, where γ-ray and neutrino emission are produced in the same location, can hardly reproduce the broadband spectral energy distribution and the temporal profiles of TXS 0506+056 simultaneously. As a result, more complicated multi-zone models are likely to be needed (Murase et al. 2018; Cerruti et al. 2019). Halzen et al. (2019) and Reimer et al. (2019) showed that the absorption and interactions intrinsic to the source, and by the extragalactic background, may significantly alter the γ-ray spectrum. Internal absorption may cause a lack of high-energy γ-ray and neutrino correlation in AGN (Aartsen et al. 2017) and maybe hint towards a correlation with lower energy observations in the X-ray and radio bands. In summary, despite multi-wavelength follow-ups of the event and extensive theoretical modeling, the radiation mechanisms of TXS 0506+056 and the underlying magnetic field strength and configuration remain poorly understood. The theory of neutrino production in blazar jets is still largely unknown.

Following the identification of TXS 0506+056 the majority of efforts have focused in the γ-ray regime. However, this connection has proven complicated. The regions that would host sufficiently dense external photon fields to produce high neutrino fluxes would also imply strong internal absorption of γ-rays (see e.g., Murase et al. 2016; Reimer et al. 2019), and therefore no temporal coincidence would be expected. Indeed, no γ-ray flare was detected from TXS 0506+056 during the 2014–2015 neutrino burst (IceCube Collaboration 2018). Recent modeling results also suggest that γ-ray emission could be heavily suppressed during neutrino flares (Kun et al. 2021; Mastichiadis & Petropoulou 2021).

Further complications in identifying neutrino counterparts arise from the possibility that not all blazar classes are equally “good” neutrino candidates (e.g., Neronov & Semikoz 2002). Padovani et al. (2016) found the most significant connection with high synchrotron peaked (HSP) sources¹. The most recent revision (with more data) suggested a connection with intermediate synchrotron peaked (ISP) and HSP sources (Giommi et al. 2020). On the other hand, Huber (2019) found the most significant connection (though only at 1.9σ level) with low synchrotron peaked (LSP) and ISP BL Lac objects, and the connection to HSP had a much lower significance (only 0.5σ). Plavin et al. (2020) carried out a careful study of 3388 jetted-AGN for which there are very long baseline interferometry (VLBI) images at 8 GHz, and a sample of 56 IceCube high-energy neutrino events E >  200 TeV, and found significant evidence for an association, with a chance probability of only 0.2%. They also found highly probable associations with the jetted-AGN 3C 279, NRAO 530, PKS 1741−038, and OR 103. In a separate study, using IceCube track data for seven years and 712 830 detected events, and 3411 jetted AGN with VLBI data, Plavin et al. (2021) found evidence for association with a post trial p-value of 3 × 10⁻³, corresponding to a 3σ significance for a normal distribution. The combined significance of these two studies is 4.1σ. Shortly after, this result was challenged by Zhou et al. (2021) using the same sample.

In a recent paper, we investigated the connection between astrophysical high-energy neutrinos and jetted-AGN monitored by the Owens valley radio observatory (OVRO) and Metsähovi radio observatory Hovatta et al. (2021, hereafter H21). OVRO has been monitoring a large sample of jetted-AGN since 2008 with regular cadence of about 3 days. The initial sample consisted of 1157 sources (hereafter CGRaBS sample). Over the years more sources detected by the Fermi gamma-ray space telescope (Fermi) were added to the sample, which at the time of H21 numbered 1795 sources (hereafter full sample). We considered three different samples: the CGRaBS sample, the full sample, and a flux-density limited sample at 350 mJy (Liodakis et al. 2017). We then cross-matched the coordinates of the OVRO sources with the reconstructed position of high-energy neutrinos, and evaluated the radio activity state of the coincident sources. Using radio light curves from OVRO, we identified 20 candidate sources positionally associated with astrophysical neutrinos during periods of enhanced radio activity at the ∼2σ level. Almost half the potentially associated sources (9/20) have not been detected by Fermi. This is interesting considering that a large fraction of the OVRO monitoring sample was included after the sources were detected by Fermi, implying that γ-ray activity does not play a major role in the neutrino association. Additionally, the fact that they have not been detected by Fermi, places strong constraints on the possible location and neutrino emission mechanism as the jet physical conditions have to be such that protons are energetic enough to produce neutrinos, but at the same time not produce any detectable γ-ray emission. Alternatively, neutrino production has to take place in regions where the opacity for γγ pair production is very high.

Apart from the 13 ± 5 neutrino excess above the atmospheric background found in TXS 0506+056, thus far all of the suggested neutrino associations only have a single high-energy neutrino detection per source (either an AGN or a tidal disruption event [TDE]). In the claimed TDE association (Stein & Velzen 2021) the neutrino lags the TDE by ∼120 days, and the association has been challenged (Cendes et al. 2021). The claimed ∼3σ AGN association has also been questioned (Böttcher 2019; Reimer et al. 2019; Kun et al. 2021). Considering different neutrino datasets yields a < 3σ probability of association with TXS 0506+056 (Abbasi et al. 2021). The a posteriori nature of these associations is not conducive to establish real connections with AGN, since a posteriori probabilities are notoriously difficult to calculate. Two of the studies based on well-defined samples that took the a priori approach showed high-energy neutrino-AGN associations significant at the 2.9σ (Plavin et al. 2020) and 3.3σ (Abbasi et al. 2021) levels.

In this paper, we focus on the potential connection between blazars and astrophysical neutrinos. We lay out the necessary strategies that will allow us to prove or disprove AGN with jets as the sources of high-energy neutrinos and explore the shortcomings of our current efforts. In Sects. 2 and 3 we lay out and discuss different strategies. In Sect. 4 we simulate neutrino counterpart samples and explore the detectability of a statistical significant correlation with jetted-AGN by current and future experiments. We draw conclusions in Sect. 5.

2. General strategies for identifying neutrino counterparts

Since we lack compelling proof that high-energy astrophysical neutrinos are associated with AGN, it is important to be as unbiased as possible in our approach. It would be just as important to disprove as to prove that there is such an association. With this in mind, we must consider how to optimize our chances. This leads us to the following conclusions regarding the optimization of any potential experiment with that goal. Firstly, it should establish clear a priori hypotheses to test, and not change the criteria once observing has begun. Clearly these hypotheses can be applied only to observations that commence after these hypotheses are laid down, that is only to data taken after the start of any experiment. In addition, sample selection should be unbiased and the sample should not try to mimic properties of existing candidates, but rather be motivated by theory. Secondly, Earth shielding against atmospheric muons is important in obtaining as clean a sample of astrophysical neutrinos as possible. Therefore this experiment should concentrate primarily on neutrinos from the hemisphere opposite to that of the neutrino detector. Thirdly, as discussed above, it is possible that not all types of blazars are good neutrino candidates, and some might be γ-ray dark. Sources that intrinsically do not produce neutrinos will dilute the sample. In this case, looking at the significance of individual sub-classes is preferable to the combined analysis. Therefore, any sample should include as many different classes of jetted-AGN as possible. Finally, the probability of finding an association with a flare is proportional to the observing time, thus at least a few years of monitoring are necessary. To the best of our knowledge only the current monitoring program by OVRO meets the above criteria.

3. Positional association alone is futile

The futility of spatial-only searches for neutrino-AGN associations becomes obvious when considering that there are ≳10³ AGN in each IceCube positional uncertainty area on the sky. We came to the same conclusion in H21. Here, we further demonstrate the need for timing information.

We consider below the statistical significance of associations using the sample of neutrinos from our work in H21 and a simple binomial calculation. Let us simplify the association of jetted-AGN and neutrinos into a “dart throwing” experiment. We have a number of N_neut neutrinos detected by IceCube, which have a certain positional uncertainty region around them. We then attempt to “hit” these neutrinos with our sample of n jetted-AGN. If the jetted-AGN falls within the uncertainty ellipse of a neutrino, we call that a success. Then the probability of getting at least k successes just by chance, that is the estimate of how unlikely k successes are, can be obtained by the binomial cumulative distribution function

$$\begin{array}{c}\hfill P(k\mathrm{or}\mathrm{more})={\displaystyle \sum _{i=k}^n}\left(\begin{array}{c}n\\ i\end{array}\right)p_{\mathrm{chance}}^i{(1-p_{\mathrm{chance}})}^{n-i},\end{array}$$(1)

where p_chance is the probability for an individual jetted-AGN to fall within the uncertainty ellipses of the N_neut neutrino events on the null hypothesis that neutrino events arise from an isotropic population.

We can estimate p_chance by using the sample of N_neut = 56 high-energy neutrinos detected by IceCube in 2009–2019 used in H21. The mean uncertainty ellipse radius, including also a systematic uncertainty of 0.9° added in quadrature, was found to be 1.5°, resulting in an ϵ = 7.1 deg² uncertainty region on the sky around each neutrino. Because OVRO only observes sources at declination δ >  −20°, we consider an area of A = 27 681 deg² on the sky.

Now the chance probability for a single jetted-AGN to fall within the uncertainty ellipses of the neutrinos is

$$\begin{array}{c}\hfill p_{\mathrm{chance}}=\frac{N_{\mathrm{neut}}ϵ}{A}=0.014.\end{array}$$(2)

To demonstrate the calculation of the binomial probability, we used the sample of 1157 CGRaBS sources from H21, which is our number of trials n. We find that 17 of these fall within the error ellipse of the N_neut = 56 neutrinos, meaning that our number of successes is k = 17. Using Eq. (1) we find that the probability to have at least this many associations by chance is P(k ≥ 17)=0.49, hence we cannot claim a connection based on positional association alone.

We now repeat the exercise for 5000 sources, assuming the same 56 neutrino events, as this is close to the expected number of events in the next 5 years. The number of sources is motivated by OVRO’s current 5-year monitoring program dedicated to the association of neutrinos with jetted-AGN. In this case, the chance probability for a single event is the same as before, p_chance = 0.014. However, now we have n = 5000 jetted-AGN, and the number of matches increases to k = 64. We note that this is actually higher than the number of observed neutrinos N_neut = 56 because multiple sources can fall within the same uncertainty ellipse. The binomial probability will be P(k ≥ 64)=0.84, which is completely uninformative.

4. Detecting a correlation between AGN and neutrinos

The above calculations clearly demonstrate how a positional association is not enough to prove the physical connection between jetted-AGN and neutrinos at a high significance. In order to do that, one needs to somehow informatively reduce the number of trials n, which can only be achieved by accounting for timing information. In H21 we showed using the CGRaBS sample that the probability for observing a strong flare from a jetted-AGN coincident in time and position with a neutrino event was p = 0.03 after accounting for all trials. This 2σ result is by no means conclusive, but certainly an improvement. Although increasing the sample size will increase the probability of chance association, it will also allow one to find more real matches, if they exist, which will improve statistics. Similar conclusions for the need of coincident flares were reached in Capel et al. (2022) looking at the neutrino – γ-ray connection through population modeling.

It is important to assess the detectability of a potential correlation with jetted-AGN given the aforementioned phenomenological challenges we face. To gain further insight in the necessary conditions and evaluate the current methodology to achieve our objectives, we simulated neutrino counterpart populations under the null hypothesis that those neutrinos are emitted during radio flares. We start by asking the simple question: can we detect a correlation if all neutrinos are astrophysical and if all were produced during the peaks of radio flares? To address this question we used the radio light curves and the neutrino sample from H21. We used the events in the H21 neutrino sample and assigned new arrival times so that neutrinos coincide with peaks of radio flares. We then repeated the same analysis as in H21. This is achieved by modeling the light curves of the associated CGRaBS sources found in H21 with Bayesian blocks (Scargle et al. 2013). We identified flares as the local maxima in the Bayesian block representation (green dash lines, Fig. 1) following Liodakis et al. (2018) and Liodakis et al. (2019). We used the width of the peak block to randomly assign an arrival time to the neutrino associated with a given source and created 1000 simulated samples. We considered three possibilities. (1) we randomly select identified flares in the light curve; (2) we randomly select flares whose peak blocks have a flux density higher than the median of the flares in that light curve; and (3) we select only the highest flux-density flare. For each of the cases, we repeated the analysis as in H21 and compared our results with the observed. The results for the different simulation cases are given in Table 1. Here we provide a short description and refer the reader to H21 for the details of the cross-correlation analysis. We first start by multiplying the reported neutrino localization RA and DEC uncertainties by 1.3 to convert them to the 90% sky containment ellipse, and then added the systematic uncertainty in quadrature. Once we have the final localization uncertainty, we cross-matched it with the coordinates of the sources in our sample to find the ones positionally associated with the neutrino region. We estimated the activity index of each “associated” source as the ratio of the mean flux-density within a 2-year interval centered around the neutrino arrival time over the mean flux-density of the entire light curve (excluding the neutrino event). We then proceeded to generate simulated samples by randomly shifting the neutrino positions only in RA to reproduce IceCube’s sensitivity declination dependence. We cross-matched the “new” neutrino positions with the observed sample and calculate the activity index for each newly associated source. We repeated this process 10⁴ times to evaluate the pre-trial chance probability of obtaining an equal or larger value than in the observed population, under the null hypothesis of no association between radio AGN and high-energy neutrinos. Finally, we estimated the post-trial (look-elsewhere effect corrected) probability by treating each simulated population as the observed and repeating the cross-matching analysis described above.

Fig. 1. Bayesian block representation (red dash-dotted line) of an observed (CGRaBS J1504+1029, top panel) and a simulated source (bottom panel). The vertical dashed green lines indicate the location of the identified flares.

The values in Cols. 4–9 of Table 1 then give the fraction of those 1000 simulations that give a post-trial p-value smaller than 3σ or 2σ for the different tests that were considered in H21. This will allow us to deduce the likelihood of obtaining a significant result under different hypotheses and conditions. We include both 3σ and 2σ limits to see when we would expect a result of higher significance than in H21 for the CGRaBS sample (3σ) and in which cases the significance is similar to what was obtained (2.3σ before accounting for the sample trial, 1.6σ if the number of samples is accounted for).

The second test considered is the number of flaring sources in the associated sample (Cols. 6 and 7). In H21, we defined a source to be in a flaring state if its activity index is larger than 1.29. In this test, we calculate how many such sources there are in the associated sample and compare this to a random population generated as explained above. The significance in H21 for the CGRaBS sample was 2.8σ before accounting for the sample trial and 2.2σ after accounting for the additional sample trial factor.

The final test (Cols. 8 and 9) calculates the mean flux density of the associated sample and compares it to the rest of the population. In H21 the significance of this test for the CGRaBS sample was 2.0σ before accounting for the sample trial and 1.2σ after accounting for the additional sample trial factor.

For the three simulated samples described above we find that if the neutrino arrives during any randomly selected flares without any flux density cutoffs (line 1 of Table 1), only a relatively small fraction of simulations would give us a 2σ result and especially in the number of flaring sources and mean flux density we would almost never obtain a 3σ result. Given that our real CGRaBS data in H21 gave us a ∼3σ result for all these tests (without considering the additional sample trials), this scenario is physically unlikely. On the other hand, if all the neutrinos were astrophysical and arrived during the peak of the highest flare in the source (line 3 of Table 1) we should have obtained a 3σ result also in H21, which was not the case so from these three tests considered, the most realistic one is the case where the neutrino arrives during a flare with an amplitude higher than the median amplitude (line 2), as this gives a relatively high fraction of 2σ results, similar to the real data in H21.

We then relaxed the assumption that all neutrinos are astrophysical. To achieve that we used the reported signalness of the events (Abbasi et al. 2022)^2,3,4. The signalness is a measure of the probability that the event is of astrophysical origins. It is defined as

$$\begin{array}{c}\hfill \mathrm{Signalness}=\frac{N_{\mathrm{S}}(E,\delta )}{N_{\mathrm{S}}(E,\delta )+N_{\mathrm{B}}(E,\delta )},\end{array}$$(3)

where N_B and N_S are the number of background and signal events at a given declination (δ) higher than a proxy-energy E (Blaufuss et al. 2019). In case a signalness value is not reported, we assumed a 50% probability of astrophysical origins. We followed the same procedure as above with the exception that when a neutrino is associated we draw a random value from a uniform distribution between [0,1]. If that value is lower than the fractional signalness we consider that neutrino to be “astrophysical” and assign it to a flare as described above. If the neutrino is deemed “atmospheric” we assigned a random arrival time within the observed light curve. Including the signalness (lines 4 and 5 of Table 1), we find that the median flare case gives a 2σ result in less than half of the cases, especially in the number of flaring sources test (Col. 7), which gave the highest significance in H21. Even in the case of the neutrino arriving during the highest flare, the fraction of 3σ results is quite small. This indicates that our sample of 1158 was insufficient to obtain a significant result in the first place.

Next we explored if a sample of 5000 sources monitored for 5 years will be sufficient to produce a correlation. The choice for the number of sources and monitoring period is again motivated by the on-going OVRO experiment⁵, which we used to build our simulated samples. We used the Timmer & Koenig method (Timmer & Koenig 1995) to produce simulated light curves following Nilsson et al. (2018). The model assumes a Gaussian probability density function of the time series and does not include flattening of the spectrum at low frequencies. We also did not include a white noise component in our model. The power spectral density (PSD) slope (β) values necessary to produce the simulated light curves were randomly drawn from the empirical cumulative distribution function taken from Max-Moerbeck et al. (2014).

About 24% of the sources in the real OVRO sample are faint (< 60 mJy) and near the limit where variability can be detected in the OVRO light curves (Richards et al. 2011). The light curves of those sources typically exhibit near white noise type variability. We assumed the PSDs of those sources follow an exponential distribution with a mean of 0.175. This results in a PSD distribution bounded between [0,1.65]. We generate ∼40 year long light curves from which we randomly select a 5-year interval. We used the observed flux densities and uncertainties for randomly selected sources to rescale our simulated light curves to the observed range. The light curves were then sampled to the 10-day cadence of the current OVRO monitoring program and the measurement uncertainties were added. An example of the final simulated light curves is shown in Fig. 1 (bottom panel).

Given the current detection rate of all “gold” neutrino events (i.e., events with > 50% probability of astrophysical origins), we expect roughly 60 events over five years (Blaufuss et al. 2019). Therefore we used the 56 neutrinos listed in H21 along with their sky positions, localization uncertainties, and signalness. We randomly assigned coordinates from the observed OVRO sample to the simulated light curves and then associated the 5000 source sample with the neutrino positions. This results in 32 neutrinos associated with the simulated sample. In 22 cases there are multiple sources within the localization uncertainty of individual neutrinos. We then repeated the analysis discussed above considering only the cases where the signalness is accounted for in the simulated samples. The results of this test are given on lines 6 and 7 of Table 1. In these cases, the fraction of 2σ results is fairly high especially when considering neutrinos arriving during the maximum flares, but it would be very unlikely for us to ever obtain a 3σ result and improve the statistics of H21 because it would also require (unrealistically) that the neutrinos arrive only during the highest flares.

So far we have assumed that neutrino emission is independent of radio flux density. To investigate a possible dependence, we replaced the “associated” sources in our simulated sample with randomly selected sources that have a median flux density of > 0.5 Jy and > 1 Jy. These tests are given on lines 8–11 of Table 1. These results show that if the neutrinos arrived during the largest flares, we would very likely get a result that is significant at the 3σ level, especially if it is the brightest radio sources that are responsible for the neutrino emission, as claimed by Plavin et al. (2020). Even the case where the neutrino arrives during brighter than median flares would give a significant result in 40% of the simulations when looking at the number of flaring sources, demonstrating that with a large enough sample it will be possible to obtain more conclusive results.

5. Conclusions

In this paper we discussed strategies that can potentially allow us to discover whether jetted-AGN are neutrino emitters. We focused on the case where all high-energy ∼PeV neutrinos are originating from jetted-AGN. In summary, monitoring of a large, statistically well defined, unbiased sample is necessary if we are to have any chance of proving or disproving a possible correlation. To further test the capabilities of our current methodology we compared the results of the observed correlation with simulated neutrino samples. We first tested a scenario with existing data used in H21 and artificially placing neutrinos during radio flares. We found that if all neutrinos are of astrophysical origin and arrive at the peak of the brightest radio flares we would have detected a significant correlation. However, we know for a fact that neutrinos do not always arrive at the peak of the flares. Placing neutrinos at any radio flare produces less significant results than in the observed case. Including the signalness in our simulations we find that even in the unreasonable case of all neutrinos coming at the peak of the brightest flares we would most likely not detect a significant correlation. This is in part due to the sample size and the number of detected/associated neutrinos, but also in part due to the uncertainty of the neutrinos’ origin. For seven out of the associated neutrinos in H21 we have assumed a 50% signalness, and three have a reported signalness of < 50%. Given that those estimates come from detector limitations, location of the conversion point etc. it is likely that they do not represent the true probability of astrophysical origin. Therefore, it is likely we are polluting our simulated samples with false negatives.

We then explored if a larger sample can produce significant results. We simulated 5000 – 5-year long – light curves which were then associated with the existing high-energy neutrino sample. We then repeated the above analysis accounting for the signalness of the neutrinos. Assuming neutrinos arrive at the peak of the brightest flares, regardless of radio brightness, then our results would not be significantly improved over H21. On the other hand, if only the brightest radio sources are neutrino emitters, then the current OVRO experiment will most likely find a significant correlation. This is particularly important given the clear relation between radio brightness and synchrotron peak (HSPs are typically < 0.5Jy). Currently, no more than one high-energy neutrino has been found within the same localization uncertainty. Detecting multiple neutrinos coming from the same sources would tremendously improve the statistical significance of our analysis. This would most likely require a longer period than the one considered for our simulations, highlighting the need for continuous long-term monitoring programs. As previously discussed, a larger sample would increase the chance coincident, but also the real associations. Unfortunately, samples > 5000 sources are beyond the monitoring capabilities of any existing single radio telescope with a reasonable cadence. Upcoming large scale radio surveys (e.g., Simons observatory – Ade et al. 2019; CMB-S4 – Abazajian et al. 2022) can provide the necessary monitoring capabilities, and would require neutrino observatories in the northern hemisphere. While our analysis is focused on the radio regime, the overall strategies are applicable to other wavelengths. Vera C. Rubin observatory’s legacy survey of space and time (LSST) will be soon monitoring millions of AGN (LSST Science Collaboration 2009). Assuming a linear relation between neutrinos and bolometric luminosity, Creque-Sarbinowski et al. (2022) found that a statistically significant correlation can be achieved for 2.8 × 10⁷ AGN monitored for the nominal 10-year duration of LSST, even if only a small fraction of AGN produce neutrinos. The authors find that the correlation significantly improves when considering AGN variability and using neutrino detectors in the opposite hemisphere, well inline with the strategies outlined here.

We note that in our analysis considering the signalness we assumed that all astrophysical neutrinos originate in jetted-AGN. If different types of sources are responsible for the high-energy neutrinos observed by IceCube, in addition to AGN (e.g., TDEs, van Velzen et al. 2021), it will have a negative impact on the significance of the correlation. From our analysis we can come to the following conclusions.

• Our simulations of neutrinos arriving at any radio flare (possibility (1) in Sect. 4) yielded a less significant correlation than that which was observed. This suggests that there is a positive relation between neutrinos and radio brightness.

• Similarly, our results imply that at least some of the neutrinos that appear to be associated with blazars are likely astrophysical.

• The sample in H21, although the largest timing study thus far, was not sufficient to yield a significant correlation. Larger samples are necessary.

• A ∼5000 source sample monitored over 5 years is sufficient to detect a correlation at a > 3σ level, only if the brightest (> 0.5 Jy) radio sources emit neutrinos.

• If all blazars are equally good neutrino emitters, it is unlikely a correlation more statistically significant than then one reported in H21 will be found.

• Although not definitive, a less significant correlation than the one reported in H21 would argue against the blazar-neutrino association.

Improving the state-of-the-art analysis tools (e.g., an improved measure for the activity index), as well as new observations from operating (ANTARES, Illuminati et al. 2021; Baikal-GVD, Belolaptikov et al. 2022) neutrino telescopes will help us mitigate the current observational limitations. Future facilities such as KM3NET (Aiello et al. 2019; Adrián-Martínez et al. 2016) as well as upgrades to IceCube (IceCube-gen2, The IceCube-Gen2 Collaboration 2021) will significantly improve the angular resolution and the accuracy of the reconstructed arrival direction. The sub-degree accuracy that could be achieved for high-energy neutrinos (> 10 TeV) will be instrumental in reducing the number of jetted-AGN within the localization uncertainty, and therefore improving the significance of spatial and temporal correlations, as well as revealing neutrino point-sources by associating multiple neutrinos to a single jetted-AGN.

Acknowledgments

The authors thank Keith Grainge and the anonymous referee for comments and suggestions that helped improve this work. T.H. was supported by the Academy of Finland projects 317383, 320085, 322535, and 345899. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), supported by private funding from the California Institute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST- 1109911. E.L. was supported by Academy of Finland projects 317636 and 320045. The computer resources of the Finnish IT Center for Science (CSC) and the FGCI project (Finland) are acknowledged. S.K. acknowledges support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 771282. W.M. gratefully acknowledges support by the ANID BASAL projects ACE210002 and FB210003, and FONDECYT 11190853. V.P. acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project 1552 CIRCE), and from the Foundation of Research and Technology – Hellas Synergy Grants Program through project MagMASim, jointly implemented by the Institute of Astrophysics and the Institute of Applied and Computational Mathematics.

References

All Tables

All Figures

	Fig. 1. Bayesian block representation (red dash-dotted line) of an observed (CGRaBS J1504+1029, top panel) and a simulated source (bottom panel). The vertical dashed green lines indicate the location of the identified flares.
In the text