Searches for Neutrinos in the Direction of Radio-bright Blazars with the ANTARES Telescope

Active galaxies, especially blazars, are among the most promising extragalactic candidates for high-energy neutrino sources. To date, ANTARES searches included these objects and used GeV–TeV γ-ray flux to select blazars. Here, a statistically complete blazar sample selected by their bright radio emission is used as the target for searches of origins of neutrinos collected by the ANTARES neutrino telescope over 13 yr of operation. The hypothesis of a neutrino–blazar directional correlation is tested by pair counting and a complementary likelihood-based approach. The resulting posttrial p-value is 3.0% (2.2σ in the two-sided convention). Additionally, a time-dependent analysis is performed to search for temporal clustering of neutrino candidates as a means of detecting neutrino flares in blazars. None of the investigated sources alone reaches a significant flare detection level. However, the presence of 18 sources with a pretrial significance above 3σ indicates a p = 1.4% (2.5σ in the two-sided convention) detection of a time-variable neutrino flux. An a posteriori investigation reveals an intriguing temporal coincidence of neutrino, radio, and γ-ray flares of the J0242+1101 blazar at a p = 0.5% (2.9σ in the two-sided convention) level. Altogether, the results presented here suggest a possible connection of neutrino candidates detected by the ANTARES telescope with radio-bright blazars.

Active galaxies, and particularly blazars, which have their jets pointed towards us, form a very promising class of neutrino source candidates [22,23,24,25].The emission from the jet is detectable across cosmological distances thanks to the relativistic beaming [26].Similar beaming effects should affect neutrinos as well, enhancing their chance to be detected on Earth.Moreover, there are hints that neutrinos are predominantly produced during major flares close to the jet origin [7,27], indicating a tight physical connection between jets and neutrino production processes.
Currently, all associations of high-energy neutrinos with celestial objects or classes thereof are fundamentally statistical.The comparison of findings based on different instruments is crucial to understand both astrophysical neutrino flux characteristics, such as its energy spectrum, and potential systematic effects inherent to their detection.Earlier searches for a neutrino signal in the ANTARES detector from the direction of blazars were performed using γ-bright blazars contained in the Fermi LAT catalog [28,29].However, recent studies suggest that emission and flares in γ-rays may not be tightly correlated with neutrino sources [30,8,31]: hadronic γrays are co-produced together with neutrinos, but quickly cascade down in energy.Synchrotron emission from blazar jets, detected on Earth as radio emission, could likely be a better tracer of relativistic beaming and activity happening close to the jet origin.
In addition, the recent publication of the IceCube Event Catalog of Alert Tracks-1 [32], containing a revised sample of neutrino candidates with a high probability of being of astrophysical origin, has been followed by a search [33] for a correlation with the population of blazars contained in the Fermi 4LAC-DR2 catalog [34] and in the Radio Fundamental Catalog1 .As no significant correlation is found, this result mitigates the previous findings of [7], and indicates that only a minority of the γ-ray bright and radio-bright AGN population can be identified as potential sources of high-energy neutrinos.
This paper presents an independent and complementary study of the possible neutrino-blazar association with data collected by the ANTARES observatory.Section 2 and Section 3 define the dataset used in the analysis: the ANTARES neutrino candidates and the radio observations of inner blazar regions.A description of the employed analysis methods and respective results follows.The hypothesis of directional correlation between neutrinos and blazar emission is tested by mean of a neutrino-blazar pair counting method (Section 4), as well as of a complementary time-integrated likelihood-based approach (Section 5).Moreover, in order to search for a time clustering of ANTARES events coming from the blazar directions, a time-dependent likelihood scan is also performed (Section 6).In Section 7, the results of a follow-up search for multimessenger time flare associations are presented.Finally, Section 8 summarizes all findings and discusses further prospects.

ANTARES detector and data sample
ANTARES [35] was an undersea high-energy neutrino telescope, located 40 km off-shore from Toulon, France, below the surface of the Mediterranean Sea.After over 15 years of operation, it took its last data in February 2022.The detector consisted of a three-dimensional array of 885 photomultiplier tubes (PMTs), distributed along 12 vertical lines of 450 m length, for a total instrumented volume of ∼ 0.01 km 3 .Each line was anchored to the seabed at a depth of about 2500 m and held taut by a buoy at the top.The PMTs collected the Cherenkov photons induced by the passage of the relativistic charged particles produced in neutrino interactions inside or near the instrumented volume.The information provided by the position, time and collected charge of the signals in the PMTs can be used to infer the direction and energy of the incident neutrino.
Two main event topologies, induced by different neutrino flavours and types of interactions, could be identified: track-like and shower-like events.Charged current (CC) interactions of muon neutrinos produce relativistic muons that can travel large distances through the medium, with Cherenkov light being induced all along the muon path, resulting in a track-like signature in the detector.The parent neutrino direction of high-quality selected tracks can be reconstructed with a median angular resolution of ∼ 0.8 • at E ν ∼ 1 TeV and below ∼ 0.4 • for E ν > 10 TeV [36] thanks to the long lever arm of this channel.Shower-like events are induced by all-flavour neutral current (NC) as well as ν e and ν τ CC interactions.This topology is characterised by an almost spherical light emission around the shower maximum, with an elongation of a few metres, which results in a worse estimation of the parent neutrino direction compared to the track channel.A median angular resolution of ∼ 3 • is achieved for high-quality selected showers with energies between 1 TeV and 0.5 PeV [36].
The data set employed in this analysis includes events recorded in ANTARES between January 29, 2007 and February 29, 2020 (3845 days of livetime), and selected using the criteria defined in [36], optimized to minimise the neutrino flux needed for a 5σ discovery of a point-like source with a ∝ E −2.0 emission spectrum.The selection includes cuts on the zenith angle, the angular error estimate and parameters describing the quality of the reconstruction.In the shower channel, an additional cut is applied on the interaction vertex, required to be located within a fiducial volume slightly larger than the instrumented volume.Recently improved calibrations have been used to reconstruct all the selected ANTARES events.This yields slightly different values of the reconstructed direction -with ∼ 98% of the selected events being reconstructed within 1 • from the previous direction -and of the quality parameters associated with each event.In addition, the estimator of the energy for track-like events has been updated, taking into account the time-evolution of the detector over the whole data-taking period.A total of 10504 track-like and 227 shower-like events survive the selection, with an expected atmospheric muon contamination of ∼ 14% and ∼ 43% for the track and the shower channel, respectively.Only the track channel is employed in the search for directional correlations, while the search for neutrino temporal flares makes use of both tracks and showers.

Radio observations of blazars
Synchrotron emission from inner regions of blazar jets is observed and measured by radio telescopes using very long baseline interferometry (VLBI) techniques.VLBI-measured flux densities represent the emission from the very central parsecs of blazars, and these measurements are used to select the objects included in this analysis.
Figure 1: Equatorial sky map showing the location of the VLBI blazars (red dots), together with the average density of ANTARES track-like and shower-like events per square degree.The event density at a given point in the sky is obtained dividing the number of events found within a cone of 10 • radius around this point by the solid angle of the cone.The surface of the markers is proportional to the VLBI 8 GHz flux density.The blazars located in the grey region at high declinations are outside the ANTARES field of view.The dashed black line shows the Galactic plane, and the Galactic center is represented as a black star marker.
A complete flux-limited sample of blazars observed by VLBI is selected using the same criteria as in [8].Accurate positions and parsec-scale flux densities of observed blazars are compiled in the Radio Fundamental Catalogue2 .Observations at 8 GHz are utilized because of their completeness across the whole sky.There are 3411 blazars with historical average flux density above 150 mJy; 3051 of them fall into the ANTARES field of view and are targeted by the searches for directional correlations.Figure 1 shows the location of the selected blazars in equatorial coordinates, together with the ANTARES neutrino candidates.As for the timedependent scan, the sources with low visibility, i.e. located above 40 • in declination, are excluded from the search, leading to a targeted catalog of 2774 blazars.The collection of observations includes geodetic VLBI [37,38,39], the Very Long Baseline Array (VLBA) calibrator surveys (VCS [40,41,42,43,44,45,46,47]), other 8 GHz global VLBI, VLBA, EVN (the European VLBI Network) and LBA (the Australian Long Baseline Array) programs [48,49,50,51,52,53,54,55].

Counting method
The first correlation analysis uses a simple method, inspired by [7], where the observable is the number of neutrino-blazar angular pairs separated by an angular distance Ψ less than x • β.Here, β is the angular uncertainty coming from the neutrino reconstruction, and x is a free parameter that varies in the interval [0.1; 2].The parameter x is meant to take into account a possible systematic difference between the output of the reconstruction algorithm and the true (unknown) angular error radius.Using Monte Carlo simulations, the relationship between the error estimate β and the true 68% containment radius Ψ 68% can be assessed.For values β ≲ 0.5, and reconstructed energies above ∼ 10 TeV, the relation between β and Ψ 68% is close to what is expected from a two-dimensional Gaussian function.At lower energies and for higher values of β, the value of Ψ 68% becomes significantly higher than the Gaussian expectation.The scan on the x parameter is then an empirical way to take into account this complicated behavior while still using the reconstruction quality information of ANTARES neutrinos on an event-by-event basis.As a consequence, the p-value obtained with this method needs a trial factor correction.
A possible correlation between the radio-flux density and neutrino emission is evaluated by performing an additional scan on the radio-flux density S 8GHz : VLBI blazars are kept in the sample if they satisfy S 8GHz > S min , and the value of S min is varied in the range [0.15; 5] Jy.
The results of the counting analysis are illustrated in Figures 2 and 3. When using the full VLBI catalog, the one-dimensional scan shows an absolute minimum for x = 0.82, where n obs = 469 pairs are observed in data, while n exp = 410.4are expected on average from random simulations (59 pairs in excess).The associated pre-trial p-value is p = 2.5 × 10 −3 (3.0σ), leading after correction to a post-trial p-value of P = 3.0 × 10 −2 (2.2σ).Note that without performing a scan, the value p(x = 1) = 5.0 × 10 −2 (2.0σ) would have been quoted instead.
The position x = 0.82 of the minimum seems to indicate that the angular uncertainty β of the neutrinos is over-estimated, which is not what is expected from Monte Carlo.However, the position of this scan minimum is subject to random fluctuations, and as the excess is not statistically significant, this observation cannot be interpreted as evidence for an issue in the ANTARES track reconstruction.
The search for a correlation with radio flux density with the two-dimensional scan (x, S min ) is shown in Figure 3.The absolute minimum is found for x = 0.82 and S min = 0.15 Jy, with a pre-trial p-value of p(S > 0.15 Jy) = 2.5 × 10 −3 , and a post-trial P (S > 0.15 Jy) = 0.26.This minimum corresponds to the previous findings of the one dimensional scan, and is obtained for the lowest value of the flux density cut, meaning that the whole VLBI blazar catalog is included.A local minimum is also visible for x = 0.42 and S min = 3.68 Jy, with a pre-trial p-value of p(S > 3.68 Jy) = 2.7 × 10 −3 .This excess is mainly driven by 3 blazars: J0609−1542 and J1743−0350, that have one very close ANTARES track (less than 0.2 • away), and J0538−4405 that has two events at 0.4 • distance.These sources are not found in the search for neutrino flares presented in Section 6, as only one neutrino is contributing for J0609−1542 and J1743−0350, and the two events located around J0538−4405 are separated by more than 6.5 years in time.
When accounting for the trial factors, the significance of this excess is P (S > 3.68 Jy) = 0.28.In addition, as can be seen in Figure 2 (dashed line), when excluding the blazars with S > 3.68 Jy from the catalog, the position of the minimum stays the same, and the pre-trial p-value is only slightly increased p = 4.5 × 10 −3 .These results indicate that the excess in the counting of neutrinos-blazar pairs is not induced by a small number of very high flux sources in the VLBI catalog.

Time-integrated likelihood analysis
A time-integrated likelihood analysis very similar to the one reported in [29] is performed, making use of more information about the ANTARES detector response than the one used in the simple counting method.The likelihood of the null hypothesis H b where only background is present is compared with an alternative hypothesis H s+b where signal events coming from blazars are present in the data.
The neutrino signal events are assumed to come from blazars in proportion to their measured VLBI flux density S 8GHz , with an energy spectrum modeled as a pure power-law E −γ .The analysis is performed for a fixed value of the spectral index γ, and repeated for values γ ∈ [1.8; 2.6] in 0.1 steps.
The likelihood is written as: where N is the total number of observed track-like events, S i is the probability density function (PDF) of the signal and B i the background one.The free parameters are the estimated number of signal µ s and background µ b events.The expression of the background-only likelihood L b is simply obtained by setting µ s = 0 in Equation 1.
The test statistics Q is defined as a likelihood ratio: where the likelihoods defined in Equation 1 are maximized with respect to the free parameters µ s and µ b .In practice, as the background-only term L b is maximal when µ b = N , the numerical procedure of likelihood maximization has to be performed only for the L s+b term.
The core ingredients of the likelihood analysis are the signal and background PDF, which are written as the product of a spatial and an energy term: where (α i , δ i ) are the equatorial coordinates, β i the angular uncertainty and E i the estimated energy of the i th neutrino candidate.The energy-dependent part of the signal PDF, g s (E), is computed using Monte Carlo simulations, by building a histogram of the reconstructed energies E of track-like events when assuming a pure power-law E −γ true distribution of the true neutrino energies E true .The background PDF, g b (E), is similarly obtained by building a histogram with Monte Carlo simulations, where the atmospheric neutrino fluxes (conventional+prompt components) and atmospheric muons fluxes are taken into account (see [56] for details).
The spatial term f b (δ i , E i ) for the background is considered to be independent of the right ascension α.Indeed, the ANTARES data set has been obtained by accumulating fifteen years of quasi-continuous measurement, therefore the non-uniformity of the detector exposure in local coordinates is averaged out by the rotation of the Earth, leading to a flat right ascension distribution of events.In practice, the value of the background spatial term f b (δ i , E i ) is estimated by a linear interpolation within a two-dimensional (sin δ, E) histogram built from Monte Carlo simulations.
The spatial signal term is obtained by summing over all the individual blazars contributions: where F(Ψ ij , δ i , β i , E i ) is the ANTARES point spread function (PSF) for track-like events, defined as the probability density for the event direction to fall within a given angular distance from the true incoming direction.The PSF is a sharply decreasing function of the space angle Ψ ij between the i th neutrino and the j th blazar, and mostly depends, by decreasing order of importance, on the angular error β, the reconstructed energy E, and the declination δ.The practical implementation of the PSF is obtained from Monte Carlo simulations, by building two-dimensional histograms in (E, β) for 8 separate bins of sin δ.
The weight of the j th blazar is proportional to its measured flux density w model j = S 8GHz , and corrected by the declination-dependent acceptance A(δ) of the ANTARES neutrino track sample (computed for the particular E −γ energy spectrum considered).For comparison, a basic scenario where all blazars have the same weight, w model j = 1, is considered.The smallest p-value obtained with the time-integrated likelihood analysis is p = 2.6 × 10 −2 (2.2σ), for a E −2.3 neutrino energy spectrum and with the radio-flux weight hypothesis.This spectral index value corresponds to the best-fit one from the ANTARES diffuse analysis [2], providing a hint that the VLBI blazars could contribute to the astrophysical neutrino diffuse flux.In comparison, when assuming an equal weight for each blazar in the likelihood, the resulting p-values ranges from p = 7.6 × 10 −2 (1.8σ) for E −2.3 up to p = 0.38 (0.9σ) for E −1.8 .The fact that smaller p-values are obtained with the flux-weight hypothesis could be an indication that the neutrino candidates correlate with blazars having a high flux density.
As the studied spatial correlation between neutrino candidates and the VLBI blazars is not statistically significant, upper limits (UL) at 90% confidence level are reported in Figure 4.The upper limits are computed for the different spectral indices of the assumed power-law neutrino flux (thin magenta lines), and their highest values as a function of the energy provides the most conservative limit curve (black dashed curve in the Figure ).For comparison, the 68% allowed values from the ANTARES diffuse analysis are plotted (with a 15% systematic on the flux normalization included as in [56]), together with results from IceCube [57].
As mentioned at the beginning of this section, our best fit for a cosmic diffuse neutrino flux (orange dotted line) has a spectral index of γ = 2.3, and can be compared to the corresponding blazar upper limit (dotted magenta line).For this particular value of γ, the ratio of the magenta and orange dotted lines is ∼ 0.2, implying that the VLBI blazars could not contribute to more than ∼ 20% of our estimated total diffuse flux of cosmic neutrinos.However, as can be seen on Figure 4, when accounting for the 68% confidence interval on the ANTARES diffuse flux estimation, the total VLBI upper-limit (black line) only weakly constrains our measurement.The For comparison, the IceCube limits [57] on Fermi 3LAC blazars for E −2.0 and E −2.2 spectra are shown in red.The thick black line shows the highest values of the single spectral index limits, and provides the most conservative upper-limit curve.For comparison, the orange-shaded band shows the 68% confidence region for the diffuse flux measured by ANTARES [2], the IceCube best fit to the astrophysical diffuse muon neutrino flux from [58] is displayed as a blue-shaded band, and the IceCube HESE spectrum from [59] is shown as blue markers.
total neutrino upper limit mildly constrains the first data point (blue marker) above 100 TeV of the IceCube HESE measurement, where the VLBI blazars would contribute ∼ 50% of the flux.

Time-dependent likelihood scan
The time-dependent scan looks for neutrino flares from the direction of the selected radio-bright blazars and also relies on an unbinned maximum likelihood method.Since this search looks for clustering of events in space and time, the knowledge of the detection time of the selected events is included in the likelihood and combined with the spatial and energy information.This is achieved by multiplying the signal and background PDFs of Equation 3 by a time-dependent term.Unlike in the time-integrated analysis, the likelihood is maximised independently at the position of each investigated source, meaning that each source is analysed separately.Therefore, the spatial signal term of Equation 4 simplifies into f s (α i , δ i , β i , E i ) = F(Ψ ij , δ i , β i , E i ).Regarding the signal time PDFs, two generic time profiles describing a temporary increase in neutrino emission -a Gaussian profile and a box profile -are tested.They are defined as: with t i being the detection time of the neutrino candidate event i, and T 0 and σ t being the unknown central time and duration of the flaring emission, respectively, both fitted in the likelihood maximisation.Given the small expected contribution of a cosmic signal in the overall data set, the background time profile is built using the time distribution of data events, ensuring a time profile proportional to the measured data.This PDF is computed applying less stringent selection criteria than those of the final sample so as to avoid statistical fluctuations, using the same approach as in [60].
At the location of each investigated source, the likelihood is maximised leaving as free parameters the number of signal events µ s , the signal spectral index γ, the central time of the flare T 0 , and the flare duration σ t .The procedure allows for a single flare per blazar to be fitted, and, for the given best-fit flare (flare with highest Q), provides the best-fit values μs , γ, T0 , σt .In the maximisation, the spectral index can take values between 1.0 and 3.5.The range includes the value predicted by the candidate acceleration mechanism (γ = 2.0) and the softer best-fit spectral indices measured by the IceCube Collaboration for the diffuse neutrino flux (γ = 2.37 [61], γ = 2.53 [62]) and for the point-like sources showing evidence for neutrino emission, TXS 0506+056 (γ = 2.1 and γ = 2.2 [4]) and NGC 1068 (γ = 3.2 [6]).Concerning the time-dependent parameters, T 0 can vary over the time range of the investigated ANTARES data (from January 29, 2007 until February 29, 2020), while σ t can take values between 1 day and 2000 days.
The test statistic Q, similarly to the time-integrated analysis, is defined as the logarithm of the ratio between the likelihood evaluated with the best-fit values of the free parameters and the background-only likelihood.However, in this case, the likelihood ratio is multiplied by a penalty term for short flares, σt ∆T , with ∆T being the allowed time range for T 0 .The purpose of the penalty term is to account for the larger trial factor that should be associated to short flares since a larger number of short flares than of long ones can be accommodated in a given time range, as described in [63].
In order to estimate the p-value of the best flare for each investigated source, the observed Q-value is compared to the test statistic distribution obtained with background-only pseudoexperiments (PEs).Each PE is obtained using sets of data randomised in right ascension to eliminate any local clustering due to potential sources keeping the corresponding source declination.In particular, the p-value is given by the fraction of background-like PEs with a value of the test statistic larger than the observed Q.The lowest obtained p-value identifies the most significant flare of the search.Finally, a trial correction that accounts for the fact that many candidates have been investigated is applied, by comparing the lowest obtained p-value to the distribution of the smallest p-values found when performing the same analysis on many background-only PEs.
Figure 5 shows the expected performance of this approach compared to that of the timeintegrated analysis, i.e. when the time information of the events is not considered.The 5σ discovery potential < n 5σ s > -defined as the mean number of injected signal events needed for a 5σ discovery in 50% of the PEs -is shown as a function of the duration of the flare.The time-dependent search performs better along almost all the investigated range of flare durations, with an improvement in the discovery potential by a factor ∼ 2 achieved for flares as short as 1 day.with the 3σ-flares found in ANTARES data (orange dots).See Table 1 for a list of these flares.The location of the Galactic plane (dashed line) and of the Galactic Center (star marker) is also shown.
The search results in 18 sources showing a flare with a pre-trial significance above 3σ for at least one of the tested time profiles.They are visualized in Figure 6 and listed in Table 1, together with the corresponding best-fit values of the free parameters.The most significant Table 1: List of radio-bright blazars for which a pre-trial significance above 3σ for at least one of the tested time profiles (Gaussian-shaped and box-shaped) has been obtained.The first three columns report the name and equatorial coordinates of the sources.The remaining columns summarise the results of the search in terms of best-fit central time of the flare T0 , flare duration σt , number of signal events μs , spectral index γ and pre-trial p-value, for the Gaussian-shaped and box-shaped signal time profile.The most significant flare found assuming each of the considered time shape is highlighted in bold.   1. Tracks (showers) are shown in blue (red).

Multi-messenger flares comparison
As a follow-up study of the findings of this analysis, the obtained best-fit neutrino flares have been compared to the radio light curves produced by the Owens Valley Radio Observatory (OVRO) [64] for those sources of Table 1 for which radio data are available.A notable overlap in time is noticed between the best-fit neutrino flare found in this analysis from the direction of J0242+1101 (PKS 0239+108) and its largest flare observed in radio, as shown in Figure 9.In view of this observation, the time distribution of the public data of the Fermi LAT γ-ray telescope and of the IceCube neutrino telescope compatible with the direction of J0242+1101 have also been studied.Figure 9 reports also the adaptive binned γ-ray light curve, obtained from Fermi data using the method described in [65,66].Remarkably, the most significant Fermi γ-ray flare for this blazar happened during the flaring emission detected in radio and the period The probability to observe 18 or more sources with a pre-trial flare significance above 3σ is 1.4%.
highlighted by the ANTARES analysis.Note that the γ-ray flare peak preceding radio is a typical scenario explained by synchrotron self-absorption of the jet base; this locates the γ-ray emission region close to the jet base [67,66].The time distribution of the IceCube track-like events in the 10-year point-source sample [68] with direction compatible with the blazar position within the 50% angular error reported by the IceCube Collaboration, is also shown.Only events with an angular uncertainty contour smaller than 10 deg 2 are depicted.While there is no evidence of time clustering of the IceCube events, a ν µ -induced track with the notable high energy of 50 TeV was detected during the flare, its reported angular uncertainty radius being 1.4 • .Furthemore, it is worth reporting that three IceCube alerts [69], namely 131165:9342044, 129933:32926212, and 128672:38561326, have a likelihood best-fit direction which lies at a distance of 1.4 • , 1.5 • , and 1.6 • , respectively, from J0242+1101.However, their detection date is subsequent to the one of the flare discussed here.The neutrino-radio-γ flare coincidence in J0242+1101 is an a-posteriori non-blind finding.It is still instructive to evaluate how likely such a correlation is to arise by chance.This study is conducted by running the whole flare fitting pipeline on background-only PEs.Each simulation of the analysis results in a fake list of neutrino flares, similar to Table 1.Then, exactly the same calculations are performed for the real list of flares and for random realizations.Specifically, blazars with the following characteristics are counted: • the global maximum of their radio/γ light curve falls within the ANTARES flare duration, T0 ± σt ; • that maximum is at least as high as for J0242+1101, compared to the median flux of the same blazar; this corresponds to maximum-to-median ratios above 1.6 for radio, and above 3.5 for γ.
Thanks to OVRO and Fermi observatories, both radio and γ-ray light curves are available for hundreds of blazars: the CGRaBS OVRO sample [27] and weekly curves from the Fermi LAT light Curve Repository3 [70].There are 335 blazars in the intersection of these samples.The outcome of this study indicates that it is relatively rare to find even a single matching blazar with the above mentioned characteristics: it appears in only p = 0.5% of random realizations; 4% for either radio or γ-rays separately.This fraction of p = 0.5% is the p-value, the chance coincidence probability to observe this kind of temporal neutrino-electromagnetic correlation.Given that the J0242+1101 analysis is a-posteriori, these findings should be considered as hints to be tested further when more data become available.

Summary
A search for associations between 3411 radio-bright blazars of a complete VLBI-flux-density limited all-sky sample, and an ANTARES neutrino event sample has been performed.
A time-integrated association search was conducted with pair-counting and likelihood-based approaches, using only the track-like neutrino candidates.The pair-counting method finds an excess of 59 neutrino-blazar pairs, with post-trial p-value of 0.03 (2.2σ).The likelihood-based method finds a similar p-value of p = 0.03, for an energy spectrum of E −2.3 and with the assumption of a neutrino flux proportional to the VLBI flux density.Upper limits on the neutrino flux from VLBI blazars are computed for different values of the spectral index E −γ , mildly constraining the IceCube HESE measurement around 100 TeV, where the VLBI blazars would contribute to ∼ 50% of the flux.
The time-dependent search for flares, conducted using an unbinned maximum likelihood method and combining track-like and shower-like samples, finds 18 sources with neutrino flares above 3σ (pre-trial).The two most significant flares come from the blazars J1355−6326 (3.7σ) and J1826+1831 (3.3σ), see Table 1.However, neither is significant alone after correcting for multiple trials.A cumulative effect is detected with p = 0.014 (2.5σ).This could be a hint of a population of blazars producing neutrinos with strong variations of the flux density.
A notable temporal correspondence between the neutrino flare from the direction of the J0242+1101 blazar, and the most prominent wide-band electromagnetic flare of that source is found.The chance probability of such a coincidence between neutrino emission, radio, and γ-ray light curves is p = 0.5%.This is a post-hoc estimate that could be a hint for a connection between neutrino and electromagnetic emission.More observational data are required to verify the correlation more robustly and look for other similar coincidences.
The first results on associating ANTARES neutrinos with radio-bright blazars are presented in this paper.Blazar neutrino emission appears highly variable and is likely to correlate with electromagnetic flares.Still, the situation is far from being clear, other classes of neutrino sources should contribute: the neutrino sky can be diverse, including a Galactic component.All neutrino observatories, IceCube, KM3NeT, and Baikal-GVD, are continuously improving to provide more detections with better directional accuracy and reliability.A strengthening of the observed correlation between blazars and neutrinos could then be expected in the coming years.Radio observational programs are ongoing to aid in reliable identification of coincidences and better understanding of associations.

Figure 2 :
Figure 2: Result of the counting analysis with the VLBI blazar catalog.The top and bottom panels show, as a function of the parameter x (defined in the text), the observed excess of pairs relative to random expectations, and the pre-trial p-value respectively.In the top panel, the blue band shows the ±1 σ confidence interval.In the bottom panel, the dashed curve shows the p-value obtained when excluding blazars with S > 3.68 Jy from the catalog.

Figure 3 :
Figure 3: Result of the two-dimensional scan over the radio-flux density S 8GHz and parameter x.The color code indicates the pre-trial p-value.

Figure 4 :
Figure 4: Upper limits on the one-flavour (ν µ + νµ ) total neutrino flux from the VLBI blazars obtained with the time-integrated likelihood analysis as a function of the neutrino energy.The ANTARES limits for each of the spectral indices tested are represented by thin solid violet lines.For comparison, the IceCube limits[57] on Fermi 3LAC blazars for E −2.0 and E −2.2 spectra are shown in red.The thick black line shows the highest values of the single spectral index limits, and provides the most conservative upper-limit curve.For comparison, the orange-shaded band shows the 68% confidence region for the diffuse flux measured by ANTARES[2], the IceCube best fit to the astrophysical diffuse muon neutrino flux from[58] is displayed as a blue-shaded band, and the IceCube HESE spectrum from[59] is shown as blue markers.

Figure 5 :
Figure 5: 5σ discovery potential in terms of mean number of signal events as a function of the simulated flare duration for the time-integrated analysis (dashed) and for the time-dependent analysis (solid).The simulated source is at a declination of δ = −40 • and the flare is centred at T 0 [MJD] = 57000.Similar results are obtained for different source declinations and central times.

Figure 6 :
Figure6: Sky map in equatorial coordinates showing the positions of the VLBI blazars coincident with the 3σ-flares found in ANTARES data (orange dots).See Table1for a list of these flares.The location of the Galactic plane (dashed line) and of the Galactic Center (star marker) is also shown.

Figure 7 :
Figure7: Weighted time distribution of the ANTARES events close to the location of J1355−6326 (top) and J1826+1831 (bottom).The dotted box and Gaussian time profiles have been drawn using the best-fit values of σt and T0 found in each case and reported in Table1.Tracks (showers) are shown in blue (red).

Figure 8 :
Figure 8: Cumulative excess of 3σ-flares found in ANTARES data.The vertical line indicates the observed count, while the histogram represents the number of sources in random realizations.The probability to observe 18 or more sources with a pre-trial flare significance above 3σ is 1.4%.

Figure 9 :
Figure 9: Multi-messenger light curves from the direction of the blazar J0242+1101 as a function of time, since 2008.Top panel: weighted time distribution of the ANTARES track-like (showerlike) events within 5 • (10 • ) from this object.The box profile has been drawn using the best-fit values of σt and T0 found in this analysis.Second panel: OVRO radio light curve.Third panel: adaptive binned γ-ray light curve obtained from Fermi LAT data.Bottom panel: weighted time distribution of the IceCube track-like events closer to J0242+1101 than their 50% angular uncertainty.The applied weight corresponds to the energy of each event.The color scale indicates the event angular distance from the source.