# Searches for neutrinos in the direction of radio-bright blazars with the ANTARES telescope

A Albert, S Alves, M André, M Ardid, S Ardid, J.-J Aubert, J Aublin, B

Baret, S Basa, Y Becherini, et al.

► To cite this version:

A Albert, S Alves, M André, M Ardid, S Ardid, et al.. Searches for neutrinos in the direction of radio-bright blazars with the ANTARES telescope. Astrophys.J., 2024, 964 (1), pp.3. 10.3847/1538-4357/ad1f5b. hal-04229327

HAL Id: hal-04229327

https://hal.science/hal-04229327

Submitted on 16 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

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

A. Albert¹,², S. Alves³, M. Andre⁴, M. Ardid⁵, S. Ardid⁵, J.-J. Aubert⁶, J. Aublin⁷, B. Baret⁷, S. Basa⁸, Y. Becherini⁷, B. Belhorma⁹, M. Bendahman⁷,10, F. Benfenati¹¹,12, V. Bertin⁶, S. Biagi¹³, M. Bissinger¹⁴, J. Boumaaza¹⁰, M. Bouta¹⁵, M. C. Bouwhuis¹⁶, H. Brânzaş¹⁷, R. Brujin¹⁶,¹⁸, J. Brunner⁶, J. Busto⁶, B. Caiffi¹⁹, D. Calvo³, S. Campion²⁰,²¹, A. Capone²⁰,²¹, L. Caramete¹⁷, F. Carenini¹¹,¹², J. Carr⁶, V. Carretero³, S. Celli²⁰,²¹, L. Cerisy⁶, M. Chabab²², R. Cherkaoui El Moursli¹⁰, T. Chiarusi¹¹, M. Circella²³, J. A. B. Coelho⁷, A. Coleiro⁷, R. Coniglione¹³, P. Coyle⁶, A. Creusot⁷, A. S. M. Cruz²⁴, A. F. Díaz²⁵, B. De Martino⁶, C. Distefano¹³, I. Di Palma²⁰,²¹, A. Domi¹⁶,¹⁸, C. Donzaud⁷,²⁶, D. Dornic⁶, D. Drouhin¹,², T. Eber¹⁴, T. van Eeden¹⁶, D. van Eijk¹⁶, S. El Hedri⁷, N. El Khayati¹⁰, A. Enzenhöfer⁶, P. Fermani²⁰,²¹, G. Ferrara¹³, F. Filippini¹¹,¹², L. Fusco²⁷, S. Gagliardini²⁰,²¹, J. García⁵, C. Gatius Oliver¹⁶, P. Gay⁷,²⁸, N. Geißelbrecht¹⁴, H. Glotin²⁹, R. Gozzini³, R. Gracia Ruiz¹⁴, K. Graf¹⁴, C. Guidi¹⁹,³⁰, L. Haegel⁷, S. Hallmann¹⁴, H. van Haren³¹, A. J. Heijboer¹⁶, Y. Hello³², J. J. Hernández-Rey³, J. Hößl¹⁴, J. Hofestädt¹⁴, F. Huang⁶, G. Illuminati¹¹, G. W. James²⁴, B. Jisse-Jung¹⁶, M. de Jong¹⁶,³³, P. de Jong¹⁶,¹⁸, M. Kadler³⁴, O. Kalekin¹⁴, U. Katz¹⁴, A. Kouchner⁷, Y. A. Kovalev³⁵,⁵⁷, Y. Y. Kovalev³⁵,³⁶,³⁷,⁵⁷, I. Kreykenbohm³⁸, V. Kulikovskiy¹⁹, R. Lahmann¹⁴, M. Lamoureux⁷, A. Lazo³, D. Lefèvre³⁹, E. Leonora⁴⁰, G. Levi¹¹,¹², S. Le Stum⁶, D. Lopez-Coto⁴¹, S. Loucatos⁷,⁴², L. Maderer⁷, J. Manczak³, M. Marcelin⁸, A. Margiotta¹¹,¹², A. Marinelli⁴³,⁴⁴, J. A. Martínez-Mora⁵, P. Migliozzi⁴⁵, A. Moussa¹⁵, R. Muller¹⁶, S. Navas⁴¹, E. Nezri⁸, B. Ó Fearraigh¹⁶, E. Oukacha⁷, A. Păun¹⁷, G. E. Păvălaş¹⁷, S. Peña-Martínez⁷, M. Perrin-Terrin⁶, V. Pestel¹⁶, P. Piattelli¹³, A. Plavin³⁵,⁴⁶,⁵⁷, C. Poiré⁵, V. Popa¹⁷, T. Pradier¹, A. Pushkarev³⁵,⁴⁷,⁵⁷, N. Randazzo⁴⁰, D. Real³, S. Reck¹⁴, G. Riccobene¹³, A. Romanov¹⁹,³⁰, A. Sánchez-Losa³,²³, A. Saina³, F. Salesa Greus³, D. F. E. Samtleben¹⁶,³³, M. Sanguineti¹⁹,³⁰, P. Sapienza¹³, J. Schnabel¹⁴, J. Schumann¹⁴, F. Schüssler⁴², J. Seneca¹⁶, M. Spurio¹¹,¹², Th. Stolarczyk⁴², M. Taiuti¹⁹,³⁰, Y. Tayalati¹⁰, S. J. Tingay²⁴, S. Troitsky⁴⁸,⁴⁹,⁵⁷, B. Vallage⁷,⁴², G. Vannoye⁶, V. Van Elewyck⁷,⁵⁰, S. Viola¹³, D. Vivolo⁴⁵,⁵¹, J. Wilms³⁸, S. Zavatarelli¹⁹, A. Zegarelli²⁰,²¹, J. D. Zornoza³, J. Zúñiga³,

(ANTARES Collaboration),

T. Hovatta⁵²,⁵³, S. Kiehlmann⁵⁴,⁵⁵, I. Liodakis⁵², V. Pavlidou⁵⁴,⁵⁵, A. C. S Readhead⁵⁶, and

(OVRO Collaboration)

¹ Université de Strasbourg, CNRS, IPHC UMR 7178, F-67000 Strasbourg, France

² Université de Haute Alsace, F-68100 Mulhouse, France

³ IFIC—Instituto de Física Corpuscular (CSIC—Universitat de València) c/ Catedrático José Beltrán, 2 E-46980 Paterna, Valencia, Spain

⁴ Technical University of Catalonia, Laboratory of Applied Bioacoustics, Rambla Exposició, 08800 Vilanova i la Geltrú, Barcelona, Spain

⁵ Institut d'Investigació per a la Gestió Integrada de les Zones Costaneres (IGIC) - Universitat Politècnica de València. C/ Paranimf 1, 46730 Gandia, Spain

⁶ Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France

⁷ Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France; julien.aublin@apc.in2p3.fr

⁸ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

⁹ National Center for Energy Sciences and Nuclear Techniques, B.P.1382, R.P.10001 Rabat, Morocco

¹⁰ University Mohammed V in Rabat, Faculty of Sciences, 4 av. Ibn Battouta, B.P. 1014, R.P. 10000, Morocco

¹¹ INFN—Sezione di Bologna, Viale Berti-Pichat 6/2, 40127 Bologna, Italy; giulia.illuminati@bo.infn.it

¹² Dipartimento di Fisica e Astronomia dell'Università di Bologna, Viale Berti-Pichat 6/2, 40127, Bologna, Italy

¹³ INFN—Laboratori Nazionali del Sud (LNS), Via S. Sofia 62, 95123 Catania, Italy

¹⁴ Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany

¹⁵ University Mohammed I, Laboratory of Physics of Matter and Radiations, B.P.717, Oujda 6000, Morocco

¹⁶ Nikhef, Science Park, Amsterdam, The Netherlands

¹⁷ Institute of Space Science, RO-077125 Bucharest, Măgurele, Romania

¹⁸ Universiteit van Amsterdam, Instituut voor Hoge-Energie Fysica, Science Park 105, 1098 XG Amsterdam, The Netherlands

¹⁹ INFN—Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

²⁰ INFN—Sezione di Roma, P.le Aldo Moro 2, 00185 Roma, Italy

²¹ Dipartimento di Fisica dell'Università La Sapienza, P.le Aldo Moro 2, 00185 Roma, Italy

²² LPHEA, Faculty of Science—Semlali, Cadi Ayyad University, P.O.B. 2390, Marrakech, Morocco.

²³ INFN—Sezione di Bari, Via E. Orabona 4, 70126 Bari, Italy

²⁴ International Centre for Radio Astronomy Research—Curtin University, Bentley, WA 6102, Australia

²⁵ Department of Computer Architecture and Technology/CITIC, University of Granada, 10807 Granada, Spain

²⁶ Université Paris-Sud, 91405 Orsay Cedex, France

²⁷ Università di Salerno e INFN Gruppo Collegato di Salerno, Dipartimento di Fisica, Via Giovanni Paolo II 132, Fisciano, 84084 Italy

²⁸ Laboratoire de Physique Corpusculaire, Clermont Université, Université Blaise Pascal, CNRS/IN2P3, BP 10448, F-63000 Clermont-Ferrand, France

²⁹ LIS, UMR Université de Toulon, Aix Marseille Université, CNRS, 83041 Toulon, France

³⁰ Dipartimento di Fisica dell'Università, Via Dodecaneso 33, 16146 Genova, Italy

³¹ Royal Netherlands Institute for Sea Research (NIOZ), Landsdiep 4, 1797 SZ 't Horntje (Texel), The Netherlands

³² Géoazur, UCA, CNRS, IRD, Observatoire de la Côte d'Azur, Sophia Antipolis, France

³³ Huygens-Kamerlingh Onnes Laboratorium, Universiteit Leiden, The Netherlands

³⁴ Institut für Theoretische Physik und Astrophysik, Universität Würzburg, Emil-Fischer Str. 31, 97074 Würzburg, Germany

³⁵ Lebedev Physical Institute of the Russian Academy of Sciences, Leninsky prospekt 53, 119991 Moscow, Russia; alexander@plav.in

³⁶ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

³⁷ Moscow Institute of Physics and Technology, Institutsky per. 9, Dolgoprudny 141700, Russia

³⁸ Dr. Remeis-Sternwarte and ECAP, Friedrich-Alexander-Universität Erlangen-Nürnberg, Sternwartstr. 7, 96049 Bamberg, Germany

³⁹ Mediterranean Institute of Oceanography (MIO), Aix-Marseille University, 13288, Marseille, Cedex 9, France; Université du Sud Toulon-Var, CNRS-INSU/IRD
UM 110, 83957, La Garde Cedex, France
# Searches for Neutrinos in the Direction of Radio-bright Blazars with the ANTARES Telescope

A. Albert<fnref n="1" />,<fnref n="2" />, S. Alves<fnref n="3" />, M. Andre<fnref n="4" />, M. Ardid<fnref n="5" />, S. Ardid<fnref n="5" />, J.-J. Aubert<fnref n="6" />, J. Aublin<fnref n="7" />, B. Baret<fnref n="7" />, S. Basa<fnref n="8" />, Y. Becherini<fnref n="7" />, B. Belhorma<fnref n="9" />, M. Bendahman<fnref n="7" />,10, F. Benfenati<fnref n="11" />,12, V. Bertin<fnref n="6" />, S. Biagi<fnref n="13" />, M. Bissinger<fnref n="14" />, J. Boumaaza<fnref n="10" />, M. Bouta<fnref n="15" />, M. C. Bouwhuis<fnref n="16" />, H. Brânzaş<fnref n="17" />, R. Brujin<fnref n="16" />,<fnref n="18" />, J. Brunner<fnref n="6" />, J. Busto<fnref n="6" />, B. Caiffi<fnref n="19" />, D. Calvo<fnref n="3" />, S. Campion<fnref n="20" />,<fnref n="21" />, A. Capone<fnref n="20" />,<fnref n="21" />, L. Caramete<fnref n="17" />, F. Carenini<fnref n="11" />,<fnref n="12" />, J. Carr<fnref n="6" />, V. Carretero<fnref n="3" />, S. Celli<fnref n="20" />,<fnref n="21" />, L. Cerisy<fnref n="6" />, M. Chabab<fnref n="22" />, R. Cherkaoui El Moursli<fnref n="10" />, T. Chiarusi<fnref n="11" />, M. Circella<fnref n="23" />, J. A. B. Coelho<fnref n="7" />, A. Coleiro<fnref n="7" />, R. Coniglione<fnref n="13" />, P. Coyle<fnref n="6" />, A. Creusot<fnref n="7" />, A. S. M. Cruz<fnref n="24" />, A. F. Díaz<fnref n="25" />, B. De Martino<fnref n="6" />, C. Distefano<fnref n="13" />, I. Di Palma<fnref n="20" />,<fnref n="21" />, A. Domi<fnref n="16" />,<fnref n="18" />, C. Donzaud<fnref n="7" />,<fnref n="26" />, D. Dornic<fnref n="6" />, D. Drouhin<fnref n="1" />,<fnref n="2" />, T. Eber<fnref n="14" />, T. van Eeden<fnref n="16" />, D. van Eijk<fnref n="16" />, S. El Hedri<fnref n="7" />, N. El Khayati<fnref n="10" />, A. Enzenhöfer<fnref n="6" />, P. Fermani<fnref n="20" />,<fnref n="21" />, G. Ferrara<fnref n="13" />, F. Filippini<fnref n="11" />,<fnref n="12" />, L. Fusco<fnref n="27" />, S. Gagliardini<fnref n="20" />,<fnref n="21" />, J. García<fnref n="5" />, C. Gatius Oliver<fnref n="16" />, P. Gay<fnref n="7" />,<fnref n="28" />, N. Geißelbrecht<fnref n="14" />, H. Glotin<fnref n="29" />, R. Gozzini<fnref n="3" />, R. Gracia Ruiz<fnref n="14" />, K. Graf<fnref n="14" />, C. Guidi<fnref n="19" />,<fnref n="30" />, L. Haegel<fnref n="7" />, S. Hallmann<fnref n="14" />, H. van Haren<fnref n="31" />, A. J. Heijboer<fnref n="16" />, Y. Hello<fnref n="32" />, J. J. Hernández-Rey<fnref n="3" />, J. Hößl<fnref n="14" />, J. Hofestädt<fnref n="14" />, F. Huang<fnref n="6" />, G. Illuminati<fnref n="11" />, G. W. James<fnref n="24" />, B. Jisse-Jung<fnref n="16" />, M. de Jong<fnref n="16" />,<fnref n="33" />, P. de Jong<fnref n="16" />,<fnref n="18" />, M. Kadler<fnref n="34" />, O. Kalekin<fnref n="14" />, U. Katz<fnref n="14" />, A. Kouchner<fnref n="7" />, Y. A. Kovalev<fnref n="35" />,⁵⁷, Y. Y. Kovalev<fnref n="35" />,<fnref n="36" />,<fnref n="37" />,⁵⁷, I. Kreykenbohm<fnref n="38" />, V. Kulikovskiy<fnref n="19" />, R. Lahmann<fnref n="14" />, M. Lamoureux<fnref n="7" />, A. Lazo<fnref n="3" />, D. Lefèvre<fnref n="39" />, E. Leonora⁴⁰, G. Levi<fnref n="11" />,<fnref n="12" />, S. Le Stum<fnref n="6" />, D. Lopez-Coto⁴¹, S. Loucatos<fnref n="7" />,⁴², L. Maderer<fnref n="7" />, J. Manczak<fnref n="3" />, M. Marcelin<fnref n="8" />, A. Margiotta<fnref n="11" />,<fnref n="12" />, A. Marinelli⁴³,⁴⁴, J. A. Martínez-Mora<fnref n="5" />, P. Migliozzi⁴⁵, A. Moussa<fnref n="15" />, R. Muller<fnref n="16" />, S. Navas⁴¹, E. Nezri<fnref n="8" />, B. Ó Fearraigh<fnref n="16" />, E. Oukacha<fnref n="7" />, A. Păun<fnref n="17" />, G. E. Păvălaş<fnref n="17" />, S. Peña-Martínez<fnref n="7" />, M. Perrin-Terrin<fnref n="6" />, V. Pestel<fnref n="16" />, P. Piattelli<fnref n="13" />, A. Plavin<fnref n="35" />,⁴⁶,⁵⁷, C. Poiré<fnref n="5" />, V. Popa<fnref n="17" />, T. Pradier<fnref n="1" />, A. Pushkarev<fnref n="35" />,⁴⁷,⁵⁷, N. Randazzo⁴⁰, D. Real<fnref n="3" />, S. Reck<fnref n="14" />, G. Riccobene<fnref n="13" />, A. Romanov<fnref n="19" />,<fnref n="30" />, A. Sánchez-Losa<fnref n="3" />,<fnref n="23" />, A. Saina<fnref n="3" />, F. Salesa Greus<fnref n="3" />, D. F. E. Samtleben<fnref n="16" />,<fnref n="33" />, M. Sanguineti<fnref n="19" />,<fnref n="30" />, P. Sapienza<fnref n="13" />, J. Schnabel<fnref n="14" />, J. Schumann<fnref n="14" />, F. Schüssler⁴², J. Seneca<fnref n="16" />, M. Spurio<fnref n="11" />,<fnref n="12" />, Th. Stolarczyk⁴², M. Taiuti<fnref n="19" />,<fnref n="30" />, Y. Tayalati<fnref n="10" />, S. J. Tingay<fnref n="24" />, S. Troitsky⁴⁸,⁴⁹,⁵⁷, B. Vallage<fnref n="7" />,⁴², G. Vannoye<fnref n="6" />, V. Van Elewyck<fnref n="7" />,⁵⁰, S. Viola<fnref n="13" />, D. Vivolo⁴⁵,⁵¹, J. Wilms<fnref n="38" />, S. Zavatarelli<fnref n="19" />, A. Zegarelli<fnref n="20" />,<fnref n="21" />, J. D. Zornoza<fnref n="3" />, J. Zúñiga<fnref n="3" />,

(ANTARES Collaboration),

T. Hovatta⁵²,⁵³, S. Kiehlmann⁵⁴,⁵⁵, I. Liodakis⁵², V. Pavlidou⁵⁴,⁵⁵, A. C. S Readhead⁵⁶, and

(OVRO Collaboration)

<footnote n="1">¹ Université de Strasbourg, CNRS, IPHC UMR 7178, F-67000 Strasbourg, France</footnote>

<footnote n="2">² Université de Haute Alsace, F-68100 Mulhouse, France</footnote>

<footnote n="3">³ IFIC—Instituto de Física Corpuscular (CSIC—Universitat de València) c/ Catedrático José Beltrán, 2 E-46980 Paterna, Valencia, Spain</footnote>

<footnote n="4">⁴ Technical University of Catalonia, Laboratory of Applied Bioacoustics, Rambla Exposició, 08800 Vilanova i la Geltrú, Barcelona, Spain</footnote>

<footnote n="5">⁵ Institut d'Investigació per a la Gestió Integrada de les Zones Costaneres (IGIC) - Universitat Politècnica de València. C/ Paranimf 1, 46730 Gandia, Spain</footnote>

<footnote n="6">⁶ Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France</footnote>

<footnote n="7">⁷ Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France; julien.aublin@apc.in2p3.fr</footnote>

<footnote n="8">⁸ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="9">⁹ National Center for Energy Sciences and Nuclear Techniques, B.P.1382, R.P.10001 Rabat, Morocco</footnote>

<footnote n="10">¹⁰ University Mohammed V in Rabat, Faculty of Sciences, 4 av. Ibn Battouta, B.P. 1014, R.P. 10000, Morocco</footnote>

<footnote n="11">¹¹ INFN—Sezione di Bologna, Viale Berti-Pichat 6/2, 40127 Bologna, Italy; giulia.illuminati@bo.infn.it</footnote>

<footnote n="12">¹² Dipartimento di Fisica e Astronomia dell'Università di Bologna, Viale Berti-Pichat 6/2, 40127, Bologna, Italy</footnote>

<footnote n="13">¹³ INFN—Laboratori Nazionali del Sud (LNS), Via S. Sofia 62, 95123 Catania, Italy</footnote>

<footnote n="14">¹⁴ Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany</footnote>

<footnote n="15">¹⁵ University Mohammed I, Laboratory of Physics of Matter and Radiations, B.P.717, Oujda 6000, Morocco</footnote>

<footnote n="16">¹⁶ Nikhef, Science Park, Amsterdam, The Netherlands</footnote>

<footnote n="17">¹⁷ Institute of Space Science, RO-077125 Bucharest, Măgurele, Romania</footnote>

<footnote n="18">¹⁸ Universiteit van Amsterdam, Instituut voor Hoge-Energie Fysica, Science Park 105, 1098 XG Amsterdam, The Netherlands</footnote>

<footnote n="19">¹⁹ INFN—Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy</footnote>

<footnote n="20">²⁰ INFN—Sezione di Roma, P.le Aldo Moro 2, 00185 Roma, Italy</footnote>

<footnote n="21">²¹ Dipartimento di Fisica dell'Università La Sapienza, P.le Aldo Moro 2, 00185 Roma, Italy</footnote>

<footnote n="22">²² LPHEA, Faculty of Science—Semlali, Cadi Ayyad University, P.O.B. 2390, Marrakech, Morocco.</footnote>

<footnote n="23">²³ INFN—Sezione di Bari, Via E. Orabona 4, 70126 Bari, Italy</footnote>

<footnote n="24">²⁴ International Centre for Radio Astronomy Research—Curtin University, Bentley, WA 6102, Australia</footnote>

<footnote n="25">²⁵ Department of Computer Architecture and Technology/CITIC, University of Granada, 10807 Granada, Spain</footnote>

<footnote n="26">²⁶ Université Paris-Sud, 91405 Orsay Cedex, France</footnote>

<footnote n="27">²⁷ Università di Salerno e INFN Gruppo Collegato di Salerno, Dipartimento di Fisica, Via Giovanni Paolo II 132, Fisciano, 84084 Italy</footnote>

<footnote n="28">²⁸ Laboratoire de Physique Corpusculaire, Clermont Université, Université Blaise Pascal, CNRS/IN2P3, BP 10448, F-63000 Clermont-Ferrand, France</footnote>

<footnote n="29">²⁹ LIS, UMR Université de Toulon, Aix Marseille Université, CNRS, 83041 Toulon, France</footnote>

<footnote n="30">³⁰ Dipartimento di Fisica dell'Università, Via Dodecaneso 33, 16146 Genova, Italy</footnote>

<footnote n="31">³¹ Royal Netherlands Institute for Sea Research (NIOZ), Landsdip 4, 1797 SZ 't Horntje (Texel), The Netherlands</footnote>

<footnote n="32">³² Géoazur, UCA, CNRS, IRD, Observatoire de la Côte d'Azur, Sophia Antipolis, France</footnote>

<footnote n="33">³³ Huygens-Kamerlingh Onnes Laboratorium, Universiteit Leiden, The Netherlands</footnote>

<footnote n="34">³⁴ Institut für Theoretische Physik und Astrophysik, Universität Würzburg, Emil-Fischer Str. 31, 97074 Würzburg, Germany</footnote>

<footnote n="35">³⁵ Lebedev Physical Institute of the Russian Academy of Sciences, Leninsky prospekt 53, 119991 Moscow, Russia; alexander@plav.in</footnote>

<footnote n="36">³⁶ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="37">³⁷ Moscow Institute of Physics and Technology, Institutsky per. 9, Dolgoprudny 141700, Russia</footnote>

<footnote n="38">³⁸ Dr. Remeis-Sternwarte and ECAP, Friedrich-Alexander-Universität Erlangen-Nürnberg, Sternwartstr. 7, 96049 Bamberg, Germany</footnote>

<footnote n="39">³⁹ Mediterranean Institute of Oceanography (MIO), Aix-Marseille University, 13288, Marseille, Cedex 9, France; Université du Sud Toulon-Var, CNRS-INSU/IRD
UM 110, 83957, La Garde Cedex, France</footnote>
⁴⁰ INFN—Sezione di Catania, Via S. Sofia 64, 95123 Catania, Italy

⁴¹ Dpto. de Física Teórica y del Cosmos & C.A.F.P.E., University of Granada, 18071 Granada, Spain

⁴² IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

⁴³ INFN—Sezione di Napoli, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy

⁴⁴ Dipartimento di Fisica Ettore Pancini, Università degli studi di Napoli Federico II, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy

⁴⁵ INFN—Sezione di Napoli, Via Cintia 80126 Napoli, Italy

⁴⁶ Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA

⁴⁷ Crimean Astrophysical Observatory, 298409 Nauchny, Crimea, Russia

⁴⁸ Institute for Nuclear Research of the Russian Academy of Sciences, 60th October Anniversary Prospect 7a, Moscow, Russia

⁴⁹ Faculty of Physics, M.V. Lomonosov Moscow State University, 1-2 Leninskie Gory, Moscow 119991, Russia

⁵⁰ Institut Universitaire de France, 75005 Paris, France

⁵¹ Dipartimento di Fisica dell'Università Federico II di Napoli, Via Cintia 80126, Napoli, Italy

⁵² Finnish Center for Astronomy with ESO, University of Turku, Vesiilinnantie 5, FI-20014, Finland

⁵³ Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland

⁵⁴ Institute of Astrophysics, Foundation for Research and Technology-Hellas, GR-70013 Heraklion, Greece

⁵⁵ Department of Physics and Institute of Theoretical and Computational Physics, University of Crete, 70013 Heraklion, Greece

⁵⁶ Owens Valley Radio Observatory, California Institute of Technology, Pasadena, CA 91125, USA

Received 2023 September 15; revised 2024 January 12; accepted 2024 January 15; published 2024 March 11

## Abstract

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 $\gamma$-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$\sigma$ 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$\sigma$ indicates a p = 1.4% (2.5$\sigma$ 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 $\gamma$-ray flares of the J0242+1101 blazar at a p = 0.5% (2.9$\sigma$ 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.

Unified Astronomy Thesaurus concepts: Neutrino astronomy (1100); Blazars (164)

## 1. Introduction

Astrophysical neutrinos of high energies, TeV and above, were discovered almost 10 yr ago by the IceCube Collaboration (Aartsen et al. 2014), consistent with later studies with the other two existing instruments, namely, ANTARES (Fusco & Versari 2019) and Baikal-GVD (Allakhverdyan et al. 2023). However, their origins are still not fully determined. Observational evidence exists for a wide range of neutrino sources, including individual active galaxies (Aartsen et al. 2018a, 2018b; Abbasi et al. 2022a), large samples of blazars (Giommi et al. 2020; Plavin et al. 2020, 2021, 2023; Buson et al. 2022, 2023; Kun et al. 2022; Abbasi et al. 2023a), and the Milky Way (Albert et al. 2017a, 2018a, 2023; Neronov & Semikoz 2016; Aartsen et al. 2019; Kovalev et al. 2022; Abbasi et al. 2023b).

Active galaxies, particularly blazars that have their jets pointed toward us, formed a very promising class of high-energy neutrino sources even before direct observational evidence became available Berezinsky (1977), Eichler (1979), Murase (2017), Boettcher (2019). The emission from the jet is

detectable across cosmological distances thanks to the relativistic beaming Blandford et al. (2019). Similar beaming effects should influence neutrinos as well, enhancing their chance to be detected on Earth (Boettcher 2023; Plavin et al. 2023). Neutrino production itself is not affected by the direction of the jet, and similar nonbeamed active galaxies can be expected to emit comparable neutrino fluxes. They are likely only detectable at closer distances, though. Moreover, there are observational hints of neutrinos being predominantly produced during major flares close to the jet origin Plavin et al. (2020), Hovatta et al. (2021), indicating a tight physical connection between relativistic 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 $\gamma$-bright blazars contained in the Fermi LAT catalog Adrian-Martinez et al. (2015); Albert et al. (2021a). However, recent studies suggest that emission and flares in $\gamma$-rays may not be tightly correlated with neutrino sources Abbasi et al. (2021), Plavin et al. (2021), Boettcher et al. (2022); hadronic $\gamma$-rays are coproduced together with neutrinos but quickly cascade down in energy. Synchrotron emission from blazar jets, detected on Earth as

<figure><img src="image_2.png" /></figure>

<footnote n="57">⁵⁷ The shown author affiliations reflect their job contracts; the ANTARES collaboration has currently suspended all institutional relations with Russian science organizations.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
⁴⁰ INFN—Sezione di Catania, Via S. Sofia 64, 95123 Catania, Italy

⁴¹ Dpto. de Física Teórica y del Cosmos & C.A.F.P.E., University of Granada, 18071 Granada, Spain

⁴² IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

⁴³ INFN—Sezione di Napoli, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy

⁴⁴ Dipartimento di Fisica Ettore Pancini, Università degli studi di Napoli Federico II, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy

⁴⁵ INFN—Sezione di Napoli, Via Cintia 80126 Napoli, Italy

⁴⁶ Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA

⁴⁷ Crimean Astrophysical Observatory, 298409 Nauchny, Crimea, Russia

⁴⁸ Institute for Nuclear Research of the Russian Academy of Sciences, 60th October Anniversary Prospect 7a, Moscow, Russia

⁴⁹ Faculty of Physics, M.V. Lomonosov Moscow State University, 1-2 Leninskie Gory, Moscow 119991, Russia

⁵⁰ Institut Universitaire de France, 75005 Paris, France

⁵¹ Dipartimento di Fisica dell'Università Federico II di Napoli, Via Cintia 80126, Napoli, Italy

⁵² Finnish Center for Astronomy with ESO, University of Turku, Vesiilinnantie 5, FI-20014, Finland

⁵³ Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland

⁵⁴ Institute of Astrophysics, Foundation for Research and Technology-Hellas, GR-70013 Heraklion, Greece

⁵⁵ Department of Physics and Institute of Theoretical and Computational Physics, University of Crete, 70013 Heraklion, Greece

⁵⁶ Owens Valley Radio Observatory, California Institute of Technology, Pasadena, CA 91125, USA

Received 2023 September 15; revised 2024 January 12; accepted 2024 January 15; published 2024 March 11

## Abstract

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 $\gamma$-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$\sigma$ 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$\sigma$ indicates a p = 1.4% (2.5$\sigma$ 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 $\gamma$-ray flares of the J0242+1101 blazar at a p = 0.5% (2.9$\sigma$ 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.

Unified Astronomy Thesaurus concepts: Neutrino astronomy (1100); Blazars (164)

## 1. Introduction

Astrophysical neutrinos of high energies, TeV and above, were discovered almost 10 yr ago by the IceCube Collaboration (Aartsen et al. 2014), consistent with later studies with the other two existing instruments, namely, ANTARES (Fusco & Versari 2019) and Baikal-GVD (Allakhverdyan et al. 2023). However, their origins are still not fully determined. Observational evidence exists for a wide range of neutrino sources, including individual active galaxies (Aartsen et al. 2018a, 2018b; Abbasi et al. 2022a), large samples of blazars (Giommi et al. 2020; Plavin et al. 2020, 2021, 2023; Buson et al. 2022, 2023; Kun et al. 2022; Abbasi et al. 2023a), and the Milky Way (Albert et al. 2017a, 2018a, 2023; Neronov & Semikoz 2016; Aartsen et al. 2019; Kovalev et al. 2022; Abbasi et al. 2023b).

Active galaxies, particularly blazars that have their jets pointed toward us, formed a very promising class of high-energy neutrino sources even before direct observational evidence became available Berezinsky (1977), Eichler (1979), Murase (2017), Boettcher (2019). The emission from the jet is

detectable across cosmological distances thanks to the relativistic beaming Blandford et al. (2019). Similar beaming effects should influence neutrinos as well, enhancing their chance to be detected on Earth (Boettcher 2023; Plavin et al. 2023). Neutrino production itself is not affected by the direction of the jet, and similar nonbeamed active galaxies can be expected to emit comparable neutrino fluxes. They are likely only detectable at closer distances, though. Moreover, there are observational hints of neutrinos being predominantly produced during major flares close to the jet origin Plavin et al. (2020), Hovatta et al. (2021), indicating a tight physical connection between relativistic 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 $\gamma$-bright blazars contained in the Fermi LAT catalog Adrian-Martinez et al. (2015); Albert et al. (2021a). However, recent studies suggest that emission and flares in $\gamma$-rays may not be tightly correlated with neutrino sources Abbasi et al. (2021), Plavin et al. (2021), Boettcher et al. (2022); hadronic $\gamma$-rays are coproduced together with neutrinos but quickly cascade down in energy. Synchrotron emission from blazar jets, detected on Earth as

<figure><img src="image_2.png" /></figure>

<footnote n="57">⁵⁷ The shown author affiliations reflect their job contracts; the ANTARES collaboration has currently suspended all institutional relations with Russian science organizations.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
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 Abbasi et al. (2023c), containing a revised sample of neutrino candidates with a high probability of being of astrophysical origin, has been followed by a search Abbasi et al. (2023d) for a correlation with the population of blazars contained in the Fermi 4LAC-DR2 catalog and the Radio Fundamental Catalog. As no significant correlation is found, this result mitigates the previous findings of Plavin et al. (2020) 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 data set 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 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.

## 2. ANTARES Detector and Data Sample

ANTARES Ageron et al. (2011) was an undersea high-energy neutrino telescope located 40 km offshore from Toulon, France, below the surface of the Mediterranean Sea. After over 15 yr of operation, it took its last data in 2022 February. 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³. 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 inside pressure-resistant optical modules Amram et al. (2002) 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 in the PMTs corresponds to a hit Aguilar et al. (2007). Hits are used to infer the direction and energy of the incident neutrino.

Two main event topologies, induced by different neutrino flavors and types of interactions, could be identified: tracklike and shower-like events Albert et al. (2021b). 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 tracklike signature in the detector Aguilar et al. (2011), Adrian-Martinez et al. (2013). The parent neutrino direction of high-quality selected tracks can be reconstructed with a median angular resolution of ~0°.8 at $E_ν \sim 1$ TeV and below ~0°.4 for $E_ν > 10$ TeV Albert et al. (2017b) thanks to the long lever arm of this channel. Shower-like events are induced by all-flavor neutral current as well as $ν_e$ and $ν_τ$ CC interactions Albert et al. (2017c). This topology is

characterized by an almost spherical light emission around the shower maximum with an elongation of a few meters, 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 Albert et al. (2017b).

The data set employed in this analysis includes events recorded in ANTARES between 2007 January 29 and 2020 February 29 (3845 days of live time) and selected using the criteria defined in Albert et al. (2017b), optimized to minimize the neutrino flux needed for a 5σ discovery of a pointlike source with a $\propto 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 the quality parameters associated with each event. In addition, the energy estimator for tracklike events has been refined to take into account the time evolution of the detector Albert et al. (2018b) over the whole data-taking period. A total of 10,504 tracklike 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.

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

A complete flux-limited sample of blazars observed by VLBI is selected using the same criteria as in Plavin et al. (2021). Accurate positions and parsec-scale flux densities of observed blazars are compiled in the Radio Fundamental Catalog.<fnref n="58" /> Observations at 8 GHz are utilized because of their completeness across the whole sky. There are 3411 blazars with a 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 time-dependent scan, the sources with low visibility, i.e., located above 40° in decl., are excluded from the search, leading to a targeted catalog of 2774 blazars. The collection of observations includes geodetic VLBI (Petrov et al. 2009; Piner et al. 2012; Pushkarev & Kovalev 2012), the Very Long Baseline Array (VLBA) calibrator surveys (Beasley et al. 2002; Fomalont et al. 2003; Petrov et al. 2005, 2006, 2008; Kovalev et al. 2007; Gordon et al. 2016; Petrov 2021), other 8 GHz global VLBI, VLBA, European VLBI Network, and Australian

<footnote n="58">58 http://astrogeo.smce.nasa.gov/sol/rfc/rfc_2020d/</footnote>
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 Abbasi et al. (2023c), containing a revised sample of neutrino candidates with a high probability of being of astrophysical origin, has been followed by a search Abbasi et al. (2023d) for a correlation with the population of blazars contained in the Fermi 4LAC-DR2 catalog and the Radio Fundamental Catalog. As no significant correlation is found, this result mitigates the previous findings of Plavin et al. (2020) 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 data set 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 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.

## 2. ANTARES Detector and Data Sample

ANTARES Ageron et al. (2011) was an undersea high-energy neutrino telescope located 40 km offshore from Toulon, France, below the surface of the Mediterranean Sea. After over 15 yr of operation, it took its last data in 2022 February. 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³. 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 inside pressure-resistant optical modules Amram et al. (2002) 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 in the PMTs corresponds to a hit Aguilar et al. (2007). Hits are used to infer the direction and energy of the incident neutrino.

Two main event topologies, induced by different neutrino flavors and types of interactions, could be identified: tracklike and shower-like events Albert et al. (2021b). 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 tracklike signature in the detector Aguilar et al. (2011), Adrian-Martinez et al. (2013). The parent neutrino direction of high-quality selected tracks can be reconstructed with a median angular resolution of ~0°.8 at $E_ν \sim 1$ TeV and below ~0°.4 for $E_ν > 10$ TeV Albert et al. (2017b) thanks to the long lever arm of this channel. Shower-like events are induced by all-flavor neutral current as well as $ν_e$ and $ν_τ$ CC interactions Albert et al. (2017c). This topology is

characterized by an almost spherical light emission around the shower maximum with an elongation of a few meters, 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 Albert et al. (2017b).

The data set employed in this analysis includes events recorded in ANTARES between 2007 January 29 and 2020 February 29 (3845 days of live time) and selected using the criteria defined in Albert et al. (2017b), optimized to minimize the neutrino flux needed for a 5σ discovery of a pointlike source with a $\propto 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 the quality parameters associated with each event. In addition, the energy estimator for tracklike events has been refined to take into account the time evolution of the detector Albert et al. (2018b) over the whole data-taking period. A total of 10,504 tracklike 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.

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

A complete flux-limited sample of blazars observed by VLBI is selected using the same criteria as in Plavin et al. (2021). Accurate positions and parsec-scale flux densities of observed blazars are compiled in the Radio Fundamental Catalog.<fnref n="58" /> Observations at 8 GHz are utilized because of their completeness across the whole sky. There are 3411 blazars with a 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 time-dependent scan, the sources with low visibility, i.e., located above 40° in decl., are excluded from the search, leading to a targeted catalog of 2774 blazars. The collection of observations includes geodetic VLBI (Petrov et al. 2009; Piner et al. 2012; Pushkarev & Kovalev 2012), the Very Long Baseline Array (VLBA) calibrator surveys (Beasley et al. 2002; Fomalont et al. 2003; Petrov et al. 2005, 2006, 2008; Kovalev et al. 2007; Gordon et al. 2016; Petrov 2021), other 8 GHz global VLBI, VLBA, European VLBI Network, and Australian

<footnote n="58">58 http://astrogeo.smce.nasa.gov/sol/rfc/rfc_2020d/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Equatorial sky map showing the location of the VLBI blazars (red dots), together with the average density of ANTARES tracklike and shower-like events per square degree. The event density at a given point in the sky is obtained by 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 gray 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 by a black star.</figcaption></figure>

Long Baseline Array programs (Petrov et al. 2011a, 2011b, 2019; Petrov 2011, 2012, 2013; Schinzel et al. 2015; Shu et al. 2017).

## 4. Counting Method

The first correlation analysis uses a simple method, inspired by Plavin et al. (2020), where the observable is the number of neutrino–blazar angular pairs separated by an angular distance $\Psi$ less than $x \cdot \beta$. Here, $\beta$ is the angular uncertainty coming from the neutrino reconstruction, and $x$ is a free parameter that varies in the interval [0.1;2.0]. 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 $\beta$ and the true 68% containment radius $\Psi_{68\%}$ can be assessed. For values $\beta \lesssim 0.5^\circ$ and reconstructed energies above $\sim 10$ TeV, the relation between $\beta$ and $\Psi_{68\%}$ is close to what is expected from a two-dimensional Gaussian function. At lower energies and for higher values of $\beta$, the value of $\Psi_{68\%}$ becomes significantly higher than the Gaussian expectation up to a factor 2, which motivated the choice of the maximum value $x=2.0$. 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_{8\text{GHz}}$; VLBI blazars are kept in the sample if they satisfy $S_{8\text{GHz}} > S_{\text{min}}$, and the value of $S_{\text{min}}$ is

varied in the range [0.15; 5.0] Jy. The maximum radio flux density of 5.0 Jy is chosen so that at least 10 sources are still included in the blazar catalog used for the correlation with neutrino events.

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_{\text{obs}} = 469$ pairs are observed in the data, while $n_{\text{exp}} = 410.4$ are expected on average from random simulations (59 pairs in excess). The associated pretrial $p$-value is $p = 2.5 \times 10^{-3}$ (3.0$\sigma$), leading after correction to a posttrial $p$-value of $P = 3.0 \times 10^{-2}$ (2.2$\sigma$). Note that without performing a scan, the value $p(x=1) = 5.0 \times 10^{-2}$ (2.0$\sigma$) would have been quoted instead.

As mentioned, the parameter $x$ is introduced to (conservatively) consider a possible systematic difference in the angular resolution not included in the reconstruction algorithm. An $x$-value smaller than 1 could represent an overestimate of the systematic uncertainty in the reconstructed track, while $x > 1$ could represent an underestimate uncertainty. The value $x = 0.82$ found for the minimum seems to indicate that the angular systematic uncertainty $\beta$ Adrian-Martinez et al. (2013) is slightly overestimated. However, as the position of this scan minimum is subject to fluctuations, this deviation of $x$ from 1 is not significant enough to conclude that the reconstruction performance behaves better than expected from Monte Carlo simulations Albert et al. (2021b).

The search for a correlation with radio flux density with the two-dimensional scan ($x, S_{\text{min}}$) is shown in Figure 3. The absolute minimum is found for $x = 0.82$ and $S_{\text{min}} = 0.15$ Jy, with a pretrial $p$-value of $p(S > 0.15 \text{ Jy}) = 2.5 \times 10^{-3}$ and a posttrial $P(S > 0.15 \text{ Jy}) = 0.26$. This minimum corresponds to
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Equatorial sky map showing the location of the VLBI blazars (red dots), together with the average density of ANTARES tracklike and shower-like events per square degree. The event density at a given point in the sky is obtained by 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 gray 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 by a black star.</figcaption></figure>

Long Baseline Array programs (Petrov et al. 2011a, 2011b, 2019; Petrov 2011, 2012, 2013; Schinzel et al. 2015; Shu et al. 2017).

## 4. Counting Method

The first correlation analysis uses a simple method, inspired by Plavin et al. (2020), where the observable is the number of neutrino–blazar angular pairs separated by an angular distance $\Psi$ less than $x \cdot \beta$. Here, $\beta$ is the angular uncertainty coming from the neutrino reconstruction, and $x$ is a free parameter that varies in the interval [0.1;2.0]. 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 $\beta$ and the true 68% containment radius $\Psi_{68\%}$ can be assessed. For values $\beta \lesssim 0.5^\circ$ and reconstructed energies above $\sim 10$ TeV, the relation between $\beta$ and $\Psi_{68\%}$ is close to what is expected from a two-dimensional Gaussian function. At lower energies and for higher values of $\beta$, the value of $\Psi_{68\%}$ becomes significantly higher than the Gaussian expectation up to a factor 2, which motivated the choice of the maximum value $x=2.0$. 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_{8\text{GHz}}$; VLBI blazars are kept in the sample if they satisfy $S_{8\text{GHz}} > S_{\text{min}}$, and the value of $S_{\text{min}}$ is

varied in the range [0.15; 5.0] Jy. The maximum radio flux density of 5.0 Jy is chosen so that at least 10 sources are still included in the blazar catalog used for the correlation with neutrino events.

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_{\text{obs}} = 469$ pairs are observed in the data, while $n_{\text{exp}} = 410.4$ are expected on average from random simulations (59 pairs in excess). The associated pretrial $p$-value is $p = 2.5 \times 10^{-3}$ (3.0$\sigma$), leading after correction to a posttrial $p$-value of $P = 3.0 \times 10^{-2}$ (2.2$\sigma$). Note that without performing a scan, the value $p(x=1) = 5.0 \times 10^{-2}$ (2.0$\sigma$) would have been quoted instead.

As mentioned, the parameter $x$ is introduced to (conservatively) consider a possible systematic difference in the angular resolution not included in the reconstruction algorithm. An $x$-value smaller than 1 could represent an overestimate of the systematic uncertainty in the reconstructed track, while $x > 1$ could represent an underestimate uncertainty. The value $x = 0.82$ found for the minimum seems to indicate that the angular systematic uncertainty $\beta$ Adrian-Martinez et al. (2013) is slightly overestimated. However, as the position of this scan minimum is subject to fluctuations, this deviation of $x$ from 1 is not significant enough to conclude that the reconstruction performance behaves better than expected from Monte Carlo simulations Albert et al. (2021b).

The search for a correlation with radio flux density with the two-dimensional scan ($x, S_{\text{min}}$) is shown in Figure 3. The absolute minimum is found for $x = 0.82$ and $S_{\text{min}} = 0.15$ Jy, with a pretrial $p$-value of $p(S > 0.15 \text{ Jy}) = 2.5 \times 10^{-3}$ and a posttrial $P(S > 0.15 \text{ Jy}) = 0.26$. This minimum corresponds to
<figure><img src="image_4.png" /><figcaption>**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 pretrial *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 (see text).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Result of the two-dimensional scan over the radio flux density $S_{8GHz}$ and parameter *x*. The color code indicates the pretrial *p*-value.</figcaption></figure>

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 pretrial *p*-value of $p(S > 3.68 \text{ Jy}) = 2.7 \times 10^{-3}$. This excess is mainly driven by three blazars: J0609−1542 and J1743−0350, which have one very close ANTARES track (less than 0.2° away), and J0538−4405, which 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 yr in time.

When accounting for the trial factors, the significance of this excess is $P(S > 3.68 \text{ 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 pretrial *p*-value is only slightly increased, $p = 4.5 \times 10^{-3}$. These results indicate that the excess in the counting of neutrino–blazar pairs is not induced by a small number of very high flux sources in the VLBI catalog.
<figure><img src="image_4.png" /><figcaption>**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 pretrial *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 (see text).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Result of the two-dimensional scan over the radio flux density $S_{8GHz}$ and parameter *x*. The color code indicates the pretrial *p*-value.</figcaption></figure>

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 pretrial *p*-value of $p(S > 3.68 \text{ Jy}) = 2.7 \times 10^{-3}$. This excess is mainly driven by three blazars: J0609−1542 and J1743−0350, which have one very close ANTARES track (less than 0.2° away), and J0538−4405, which 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 yr in time.

When accounting for the trial factors, the significance of this excess is $P(S > 3.68 \text{ 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 pretrial *p*-value is only slightly increased, $p = 4.5 \times 10^{-3}$. These results indicate that the excess in the counting of neutrino–blazar pairs is not induced by a small number of very high flux sources in the VLBI catalog.
## 5. Time-integrated Likelihood Analysis

A time-integrated likelihood analysis very similar to the one reported in Albert et al. (2021a) 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^{-\gamma}$. The analysis is performed for a fixed value of the spectral index $\gamma$ and repeated for values $\gamma \in [1.8; 2.6]$ in 0.1 steps. This range includes the value predicted by the candidate acceleration mechanism ($\gamma = 2.0$) and the softer best-fit spectral indices measured by the IceCube Collaboration for the diffuse neutrino flux ($\gamma = 2.37$, Abbasi et al. 2022b; $\gamma = 2.53$, Aartsen et al. 2020) and for TXS 0506+056, a blazar showing evidence for neutrino emission ($\gamma = 2.1$ and $2.2$; Aartsen et al. 2018a).

The likelihood is written as

$$ \ln \mathcal{L}_{s+b} = \sum_{i}^{N} \ln (\mu_s S_i + \mu_b B_i) - \mu_s - \mu_b, \quad (1) $$

where $N$ is the total number of observed tracklike events, $S_i$ is the probability density function (PDF) of the signal, and $B_i$ is the background one. The free parameters are the estimated number of signal, $\mu_s$, and background, $\mu_b$, events. The expression of the background-only likelihood $\mathcal{L}_b$ is simply obtained by setting $\mu_s = 0$ in Equation (1).

The test statistics $Q$ is defined as a likelihood ratio,

$$ Q = 2 \ln \left[ \frac{\max(\mathcal{L}_{s+b})}{\max(\mathcal{L}_b)} \right], \quad (2) $$

where the likelihoods defined in Equation (1) are maximized with respect to the free parameters $\mu_s$ and $\mu_b$. In practice, as the background-only term $\mathcal{L}_b$ is maximal when $\mu_b = N$, the numerical procedure of likelihood maximization has to be performed only for the $\mathcal{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,

$$ S_i = f_s(\alpha_i, \delta_i, \beta_i, E_i) g_s(E_i) \quad \text{and} \quad B_i = f_b(\delta_i, E_i) g_b(E_i), \quad (3) $$

where $(\alpha_i, \delta_i)$ are the equatorial coordinates, $\beta_i$ is the angular uncertainty, and $E_i$ is 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 tracklike events when assuming a pure power-law $E^{-\gamma}$ 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 muon fluxes are taken into account (see Albert et al. 2018c for details).

The spatial term $f_b(\delta_i, E_i)$ for the background is considered to be independent of the R.A. $\alpha$. Indeed, the ANTARES data set has been obtained by accumulating 15 yr of quasi-continuous measurement; therefore, the nonuniformity of the detector

exposure in local coordinates is averaged out by the rotation of the Earth, leading to a flat R.A. distribution of events. In practice, the value of the background spatial term $f_b(\delta_i, E_i)$ is estimated by a linear interpolation within a two-dimensional ($\sin \delta, E$) histogram built from Monte Carlo simulations.

The spatial signal term is obtained by summing over all the individual blazar contributions,

$$ f_s(\alpha_i, \delta_i, \beta_i, E_i) = \frac{1}{\sum w_j} \sum_{j=1}^{N_{sources}} w_j \mathcal{F}(\Psi_{ij}, \delta_i, \beta_i, E_i) \quad w_j = w_j^{\text{model}} \mathcal{A}(\delta_j), \quad (4) $$

where $\mathcal{F}(\Psi_{ij}, \delta_i, \beta_i, E_i)$ is the ANTARES point-spread function (PSF) for tracklike 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 $\Psi_{ij}$ between the $i$th neutrino and the $j$th blazar and mostly depends, by decreasing order of importance, on the angular error $\beta$, the reconstructed energy $E$, and the decl. $\delta$. The practical implementation of the PSF is obtained from Monte Carlo simulations by building two-dimensional histograms in $(E, \beta)$ for eight separate bins of $\sin \delta$.

The weight of the $j$th blazar is proportional to its measured flux density $w_j^{\text{model}} = S_{8GHz}$ and corrected by the decl.-dependent acceptance $\mathcal{A}(\delta)$ of the ANTARES neutrino track sample (computed for the particular $E^{-\gamma}$ energy spectrum considered). For comparison, a basic scenario where all blazars have the same weight, $w_j^{\text{model}} = 1$, is considered.

The smallest $p$-value obtained with the time-integrated likelihood analysis is $p = 2.6 \times 10^{-2}$ ($2.2\sigma$) for the radio flux weighting hypothesis and $p = 7.6 \times 10^{-2}$ ($1.8\sigma$) for the equal weighting. These best fits are obtained with an $E^{-2.3}$ neutrino energy spectrum for both weighting schemes. This spectral index is also the best-fit value obtained with the latest ANTARES diffuse analysis Fusco & Versari (2019). This indicates that the candidate neutrinos in the direction of the VLBI blazars do not show a different energy spectrum than the diffuse astrophysical component extracted from the observed neutrino spectrum in the ANTARES data set.

As the studied spatial correlation between neutrino candidates and the VLBI blazars is not highly significant, upper limits at a 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 violet lines), and their highest values as a function of the energy provide the most conservative limit curve (black solid 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 Albert et al. 2018c), together with the latest diffuse results from IceCube Abbasi et al. (2022c).

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 $\gamma = 2.3$ and can be compared to the corresponding blazar upper limit (black line). For this particular value of $\gamma$, the ratio between the black and orange dotted lines is $\sim 0.2$, implying that the VLBI blazars could not contribute to more than $\sim 20\%$ of our estimated total diffuse flux of cosmic neutrinos. However, as can be seen in Figure 4, when accounting for the wide 68% confidence interval on the ANTARES diffuse flux estimation, the total VLBI upper limit
## 5. Time-integrated Likelihood Analysis

A time-integrated likelihood analysis very similar to the one reported in Albert et al. (2021a) 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^{-\gamma}$. The analysis is performed for a fixed value of the spectral index $\gamma$ and repeated for values $\gamma \in [1.8; 2.6]$ in 0.1 steps. This range includes the value predicted by the candidate acceleration mechanism ($\gamma = 2.0$) and the softer best-fit spectral indices measured by the IceCube Collaboration for the diffuse neutrino flux ($\gamma = 2.37$, Abbasi et al. 2022b; $\gamma = 2.53$, Aartsen et al. 2020) and for TXS 0506+056, a blazar showing evidence for neutrino emission ($\gamma = 2.1$ and $2.2$; Aartsen et al. 2018a).

The likelihood is written as

$$ \ln \mathcal{L}_{s+b} = \sum_{i}^{N} \ln (\mu_s S_i + \mu_b B_i) - \mu_s - \mu_b, \quad (1) $$

where $N$ is the total number of observed tracklike events, $S_i$ is the probability density function (PDF) of the signal, and $B_i$ is the background one. The free parameters are the estimated number of signal, $\mu_s$, and background, $\mu_b$, events. The expression of the background-only likelihood $\mathcal{L}_b$ is simply obtained by setting $\mu_s = 0$ in Equation (1).

The test statistics $Q$ is defined as a likelihood ratio,

$$ Q = 2 \ln \left[ \frac{\max(\mathcal{L}_{s+b})}{\max(\mathcal{L}_b)} \right], \quad (2) $$

where the likelihoods defined in Equation (1) are maximized with respect to the free parameters $\mu_s$ and $\mu_b$. In practice, as the background-only term $\mathcal{L}_b$ is maximal when $\mu_b = N$, the numerical procedure of likelihood maximization has to be performed only for the $\mathcal{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,

$$ S_i = f_s(\alpha_i, \delta_i, \beta_i, E_i) g_s(E_i) \quad \text{and} \quad B_i = f_b(\delta_i, E_i) g_b(E_i), \quad (3) $$

where $(\alpha_i, \delta_i)$ are the equatorial coordinates, $\beta_i$ is the angular uncertainty, and $E_i$ is 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 tracklike events when assuming a pure power-law $E^{-\gamma}$ 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 muon fluxes are taken into account (see Albert et al. 2018c for details).

The spatial term $f_b(\delta_i, E_i)$ for the background is considered to be independent of the R.A. $\alpha$. Indeed, the ANTARES data set has been obtained by accumulating 15 yr of quasi-continuous measurement; therefore, the nonuniformity of the detector

exposure in local coordinates is averaged out by the rotation of the Earth, leading to a flat R.A. distribution of events. In practice, the value of the background spatial term $f_b(\delta_i, E_i)$ is estimated by a linear interpolation within a two-dimensional ($\sin \delta, E$) histogram built from Monte Carlo simulations.

The spatial signal term is obtained by summing over all the individual blazar contributions,

$$ f_s(\alpha_i, \delta_i, \beta_i, E_i) = \frac{1}{\sum w_j} \sum_{j=1}^{N_{sources}} w_j \mathcal{F}(\Psi_{ij}, \delta_i, \beta_i, E_i) \quad w_j = w_j^{\text{model}} \mathcal{A}(\delta_j), \quad (4) $$

where $\mathcal{F}(\Psi_{ij}, \delta_i, \beta_i, E_i)$ is the ANTARES point-spread function (PSF) for tracklike 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 $\Psi_{ij}$ between the $i$th neutrino and the $j$th blazar and mostly depends, by decreasing order of importance, on the angular error $\beta$, the reconstructed energy $E$, and the decl. $\delta$. The practical implementation of the PSF is obtained from Monte Carlo simulations by building two-dimensional histograms in $(E, \beta)$ for eight separate bins of $\sin \delta$.

The weight of the $j$th blazar is proportional to its measured flux density $w_j^{\text{model}} = S_{8GHz}$ and corrected by the decl.-dependent acceptance $\mathcal{A}(\delta)$ of the ANTARES neutrino track sample (computed for the particular $E^{-\gamma}$ energy spectrum considered). For comparison, a basic scenario where all blazars have the same weight, $w_j^{\text{model}} = 1$, is considered.

The smallest $p$-value obtained with the time-integrated likelihood analysis is $p = 2.6 \times 10^{-2}$ ($2.2\sigma$) for the radio flux weighting hypothesis and $p = 7.6 \times 10^{-2}$ ($1.8\sigma$) for the equal weighting. These best fits are obtained with an $E^{-2.3}$ neutrino energy spectrum for both weighting schemes. This spectral index is also the best-fit value obtained with the latest ANTARES diffuse analysis Fusco & Versari (2019). This indicates that the candidate neutrinos in the direction of the VLBI blazars do not show a different energy spectrum than the diffuse astrophysical component extracted from the observed neutrino spectrum in the ANTARES data set.

As the studied spatial correlation between neutrino candidates and the VLBI blazars is not highly significant, upper limits at a 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 violet lines), and their highest values as a function of the energy provide the most conservative limit curve (black solid 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 Albert et al. 2018c), together with the latest diffuse results from IceCube Abbasi et al. (2022c).

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 $\gamma = 2.3$ and can be compared to the corresponding blazar upper limit (black line). For this particular value of $\gamma$, the ratio between the black and orange dotted lines is $\sim 0.2$, implying that the VLBI blazars could not contribute to more than $\sim 20\%$ of our estimated total diffuse flux of cosmic neutrinos. However, as can be seen in Figure 4, when accounting for the wide 68% confidence interval on the ANTARES diffuse flux estimation, the total VLBI upper limit
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Upper limits on the one-flavor ($\nu_{\mu} + \bar{\nu}_{\mu}$) total neutrino flux from VLBI blazars, obtained through time-integrated likelihood analysis as a function of neutrino energy. Thin solid violet lines represent ANTARES limits for each tested spectral index (the thick black line highlights the $E^{-2.3}$ case), while the thick violet line depicts the uppermost envelope of individual spectral index limits, providing the most conservative upper-limit curve. For comparison, upper limits using IceCube data are shown (Zhou et al. 2021, thick red line; Abbasi et al. 2023d, thick blue line), together with the lower limit estimated in Plavin et al. (2021). The orange-shaded band represents the 68% confidence region for the diffuse flux measured by ANTARES (Fusco & Versari 2019), and the most recent IceCube diffuse ($\nu_{\mu} + \bar{\nu}_{\mu}$) flux estimation from Abbasi et al. (2022c) is displayed as a blue shaded band for the single power-law $E^{-2.37}$ hypothesis.</figcaption></figure>

(violet line) only weakly constrains our measurement. When comparing with the most recent IceCube diffuse ($\nu_{\mu} + \bar{\nu}_{\mu}$) flux estimation from Abbasi et al. (2022c), our most conservative limit represents ~65% of the integrated energy flux measured by IceCube between 15 and 315 TeV (the energy range where there is an overlap).

To put our constraints of blazars' contribution to the neutrino flux into context, estimates from analyses based on IceCube measurements are also shown in Figure 4. These include upper limits obtained with different versions of the VLBI catalog, rfc2019c for Zhou et al. (2021) and rfc2022a for Abbasi et al. (2023d), and lower limits estimated in Plavin et al. (2021) with version rfc2019c. The apparent tension between these estimates stems from differences in the assumptions put into corresponding statistical analyses; see discussions in Abbasi et al. (2023d) and Plavin et al. (2023). Utilizing observational data from other neutrino telescopes is crucial for independent checks of such assumptions, especially considering the statistical nature of all current associations.

## 6. 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 maximized independently at the position of each investigated source, meaning that each source is analyzed separately. Therefore, the spatial signal term of Equation (4) simplifies into $f_s(\alpha_i, \delta_i, \beta_i, E_i) = \mathcal{F}(\Psi_{ij}, \delta_i, \beta_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

$$ 
\begin{aligned}
S_{\text{Gaussian}}(t_i) &= \frac{1}{\sqrt{2\pi}\sigma_t} e^{-\frac{(t_i-T_0)^2}{2\sigma_t^2}}, & S_{\text{box}}(t_i) &= \\
&= \begin{cases} \frac{1}{2\sigma_t}, & \text{if } [T_0 - \sigma_t] \le t_i \le [T_0 + \sigma_t]; \\ 0, & \text{otherwise;} \end{cases} & (5)
\end{aligned}
$$

with $t_i$ being the detection time of the neutrino candidate event $i$, and $T_0$ and $\sigma_t$ being the unknown central time and duration of the flaring emission, respectively, both fitted in the likelihood maximization. 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 by applying less stringent selection criteria than those of the final sample so as to avoid statistical fluctuations, using the same approach as in Albert et al. (2019).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Upper limits on the one-flavor ($\nu_{\mu} + \bar{\nu}_{\mu}$) total neutrino flux from VLBI blazars, obtained through time-integrated likelihood analysis as a function of neutrino energy. Thin solid violet lines represent ANTARES limits for each tested spectral index (the thick black line highlights the $E^{-2.3}$ case), while the thick violet line depicts the uppermost envelope of individual spectral index limits, providing the most conservative upper-limit curve. For comparison, upper limits using IceCube data are shown (Zhou et al. 2021, thick red line; Abbasi et al. 2023d, thick blue line), together with the lower limit estimated in Plavin et al. (2021). The orange-shaded band represents the 68% confidence region for the diffuse flux measured by ANTARES (Fusco & Versari 2019), and the most recent IceCube diffuse ($\nu_{\mu} + \bar{\nu}_{\mu}$) flux estimation from Abbasi et al. (2022c) is displayed as a blue shaded band for the single power-law $E^{-2.37}$ hypothesis.</figcaption></figure>

(violet line) only weakly constrains our measurement. When comparing with the most recent IceCube diffuse ($\nu_{\mu} + \bar{\nu}_{\mu}$) flux estimation from Abbasi et al. (2022c), our most conservative limit represents ~65% of the integrated energy flux measured by IceCube between 15 and 315 TeV (the energy range where there is an overlap).

To put our constraints of blazars' contribution to the neutrino flux into context, estimates from analyses based on IceCube measurements are also shown in Figure 4. These include upper limits obtained with different versions of the VLBI catalog, rfc2019c for Zhou et al. (2021) and rfc2022a for Abbasi et al. (2023d), and lower limits estimated in Plavin et al. (2021) with version rfc2019c. The apparent tension between these estimates stems from differences in the assumptions put into corresponding statistical analyses; see discussions in Abbasi et al. (2023d) and Plavin et al. (2023). Utilizing observational data from other neutrino telescopes is crucial for independent checks of such assumptions, especially considering the statistical nature of all current associations.

## 6. 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 maximized independently at the position of each investigated source, meaning that each source is analyzed separately. Therefore, the spatial signal term of Equation (4) simplifies into $f_s(\alpha_i, \delta_i, \beta_i, E_i) = \mathcal{F}(\Psi_{ij}, \delta_i, \beta_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

$$ 
\begin{aligned}
S_{\text{Gaussian}}(t_i) &= \frac{1}{\sqrt{2\pi}\sigma_t} e^{-\frac{(t_i - T_0)^2}{2\sigma_t^2}}, & S_{\text{box}}(t_i) &= \\
&= \begin{cases} \frac{1}{2\sigma_t}, & \text{if } [T_0 - \sigma_t] \le t_i \le [T_0 + \sigma_t]; \\ 0, & \text{otherwise;} \end{cases} & (5)
\end{aligned}
$$

with $t_i$ being the detection time of the neutrino candidate event $i$, and $T_0$ and $\sigma_t$ being the unknown central time and duration of the flaring emission, respectively, both fitted in the likelihood maximization. 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 by applying less stringent selection criteria than those of the final sample so as to avoid statistical fluctuations, using the same approach as in Albert et al. (2019).
<figure><img src="image_7.png" /><figcaption>**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 decl. of $δ = -40°$, and the flare is centered at $T_0$ [MJD] = 57000. Similar results are obtained for different source declinations and central times.</figcaption></figure>

At the location of each investigated source, the likelihood is maximized, leaving as free parameters the number of signal events $\mu_s$, the signal spectral index $\gamma$, the central time of the flare $T_0$, and the flare duration $\sigma_t$. The procedure allows for a single flare per blazar to be fitted and, for the given best-fit flare (the flare with the highest $Q$), provides the best-fit values $\hat{\mu}_s$, $\hat{\gamma}$, $\hat{T}_0$, and $\hat{\sigma}_t$. In the maximization, the spectral index can take values between 1.0 and 3.5. Concerning the time-dependent parameters, $T_0$ can vary over the time range of the investigated ANTARES data (from 2007 January 29 to 2020 February 29), while $\sigma_t$ can take values between 1 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, $\hat{\sigma}_t/\Delta T$, with $\Delta 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 with short flares since a larger number of short flares than of long ones can be accommodated in a given time range, as described in Braun et al. (2010).

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 pseudo-experiments (PEs). In each PE, the analysis is performed on the data set randomized in R.A. to eliminate any local clustering due to potential sources keeping the corresponding source decl. 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 time-integrated analysis, i.e., when the time information of the events is not considered. The 5σ discovery potential $\langle n_s^{5\sigma} \rangle$—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 the entire investigated range of flare durations, with an improvement in the discovery potential by a factor of ~2 achieved for flares as short as 1 day.

The search results in 18 sources showing a flare with a pretrial 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 Gaussian (box) flare is found from the direction of J1355−6326 (correspondingly, J1826+1831), with a pretrial significance of 3.7σ (correspondingly, 3.3σ) in the two-sided convention. When accounting for the fact that multiple sources have been investigated, posttrial *p*-values of 29% and 84% are obtained for the best Gaussian and box flares, respectively. The result found for J1355−6326, with a best-fit flare duration of over 1000 days, is expected to be very close to the one that would have been obtained in a search for neutrino clustering in which the event time is not considered, as it can be derived from Figure 5. The weighted time distribution of the ANTARES events close to J1355−6326 and J1826+1831 is shown in Figure 7. Only tracklike (shower-like) events within a distance of 5° (10°) from the blazars are included in the plot. A higher weight is associated with events with smaller distances to the source and larger values of the energy estimator. Since the weight does not take into account the arrival time of the events, the best-fit flares do not necessarily cover in time the individual neutrino candidates
<figure><img src="image_7.png" /><figcaption>**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 decl. of $δ = -40°$, and the flare is centered at $T_0$ [MJD] = 57000. Similar results are obtained for different source declinations and central times.</figcaption></figure>

At the location of each investigated source, the likelihood is maximized, leaving as free parameters the number of signal events $\mu_s$, the signal spectral index $\gamma$, the central time of the flare $T_0$, and the flare duration $\sigma_t$. The procedure allows for a single flare per blazar to be fitted and, for the given best-fit flare (the flare with the highest $Q$), provides the best-fit values $\hat{\mu}_s$, $\hat{\gamma}$, $\hat{T}_0$, and $\hat{\sigma}_t$. In the maximization, the spectral index can take values between 1.0 and 3.5. Concerning the time-dependent parameters, $T_0$ can vary over the time range of the investigated ANTARES data (from 2007 January 29 to 2020 February 29), while $\sigma_t$ can take values between 1 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, $\hat{\sigma}_t/\Delta T$, with $\Delta 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 with short flares since a larger number of short flares than of long ones can be accommodated in a given time range, as described in Braun et al. (2010).

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 pseudo-experiments (PEs). In each PE, the analysis is performed on the data set randomized in R.A. to eliminate any local clustering due to potential sources keeping the corresponding source decl. 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 time-integrated analysis, i.e., when the time information of the events is not considered. The 5σ discovery potential $\langle n_s^{5\sigma} \rangle$—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 the entire investigated range of flare durations, with an improvement in the discovery potential by a factor of ~2 achieved for flares as short as 1 day.

The search results in 18 sources showing a flare with a pretrial significance above 3$\sigma$ 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 Gaussian (box) flare is found from the direction of J1355−6326 (correspondingly, J1826+1831), with a pretrial significance of 3.7$\sigma$ (correspondingly, 3.3$\sigma$) in the two-sided convention. When accounting for the fact that multiple sources have been investigated, posttrial $p$-values of 29% and 84% are obtained for the best Gaussian and box flares, respectively. The result found for J1355−6326, with a best-fit flare duration of over 1000 days, is expected to be very close to the one that would have been obtained in a search for neutrino clustering in which the event time is not considered, as it can be derived from Figure 5. The weighted time distribution of the ANTARES events close to J1355−6326 and J1826+1831 is shown in Figure 7. Only tracklike (shower-like) events within a distance of 5° (10°) from the blazars are included in the plot. A higher weight is associated with events with smaller distances to the source and larger values of the energy estimator. Since the weight does not take into account the arrival time of the events, the best-fit flares do not necessarily cover in time the individual neutrino candidates
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Sky map in equatorial coordinates showing the positions of the VLBI blazars coincident with the $3\sigma$ flares found in ANTARES data (orange dots). See Table 1 for a list of these flares. The two blazars with the most significant excesses are depicted in orange, while the red dot indicates the blazar studied in Section 7. The locations of the Galactic plane (dashed line) and Galactic center (star) are also shown.</figcaption></figure>

<table><caption>**Table 1**
Results of the Time-dependent Likelihood Scan</caption><thead><tr><th colspan="3">Source</th><th colspan="6">Results</th></tr><tr><th rowspan="2">Name</th><th colspan="2">δ</th><th colspan="5">Gaussian-shaped Time Profile</th><th colspan="5">Box-shaped Time Profile</th></tr><tr><th>(deg)</th><th>α</th><th>τ̂<sub>0</sub><br>(MJD)</th><th>δ̂<sub>t</sub><br>(days)</th><th>μ̂<sub>s</sub></th><th>γ̂</th><th>p-value</th><th>τ̂<sub>0</sub><br>(MJD)</th><th>δ̂<sub>t</sub><br>(days)</th><th>μ̂<sub>s</sub></th><th>γ̂</th><th>p-value</th></tr></thead><tbody><tr><td>J0112-6634</td><td>-66.6</td><td>18.1</td><td>58215</td><td>304</td><td>4.5</td><td>2.7</td><td>0.0026</td><td>58154</td><td>305</td><td>3.8</td><td>2.7</td><td>0.0097</td></tr><tr><td><b>J1355-6326</b></td><td><b>-63.4</b></td><td><b>208.9</b></td><td><b>56524</b></td><td><b>1041</b></td><td><b>7.9</b></td><td><b>2.8</b></td><td><b>0.00018</b></td><td>56091</td><td>905</td><td>6.0</td><td>2.9</td><td>0.0048</td></tr><tr><td>J0359-6154</td><td>-61.9</td><td>59.7</td><td>56316</td><td>78</td><td>5.4</td><td>3.5</td><td>0.0022</td><td>56321</td><td>112</td><td>5.7</td><td>3.5</td><td>0.0013</td></tr><tr><td>J0522-6107</td><td>-61.1</td><td>80.6</td><td>56221</td><td>42</td><td>4.6</td><td>3.4</td><td>0.0034</td><td>56232</td><td>59</td><td>4.9</td><td>3.4</td><td>0.0023</td></tr><tr><td>J1220-5604</td><td>-56.1</td><td>185.1</td><td>58406</td><td>18</td><td>2.8</td><td>2.6</td><td>0.00029</td><td>58413</td><td>22</td><td>0.4</td><td>2.2</td><td>0.0032</td></tr><tr><td>J1825-5230</td><td>-52.5</td><td>276.3</td><td>57265</td><td>600</td><td>5.4</td><td>2.7</td><td>0.0031</td><td>57188</td><td>959</td><td>5.5</td><td>2.7</td><td>0.0027</td></tr><tr><td>J0641-3554</td><td>-35.9</td><td>100.3</td><td>58084</td><td>16</td><td>2.9</td><td>3.0</td><td>0.0021</td><td>58081</td><td>19</td><td>3.0</td><td>3.0</td><td>0.0018</td></tr><tr><td>J1418-3509</td><td>-35.2</td><td>214.7</td><td>58120</td><td>11</td><td>3.3</td><td>2.9</td><td>0.0018</td><td>58121</td><td>14</td><td>2.9</td><td>2.8</td><td>0.0021</td></tr><tr><td>J1500-2358</td><td>-24.0</td><td>225.2</td><td>55846</td><td>4</td><td>3.7</td><td>2.3</td><td>0.0016</td><td>55847</td><td>6</td><td>3.7</td><td>2.2</td><td>0.0015</td></tr><tr><td>J0521-1737</td><td>-17.6</td><td>80.3</td><td>57332</td><td>1</td><td>2.0</td><td>1.9</td><td>0.0011</td><td>57333</td><td>1</td><td>2.0</td><td>1.9</td><td>0.0023</td></tr><tr><td>J2345-1555</td><td>-15.9</td><td>356.3</td><td>57653</td><td>460</td><td>3.2</td><td>2.6</td><td>0.0011</td><td>57784</td><td>404</td><td>2.4</td><td>2.7</td><td>0.0030</td></tr><tr><td>J1537-1259</td><td>-13.0</td><td>234.3</td><td>58201</td><td>46</td><td>2.6</td><td>2.0</td><td>0.0019</td><td>58201</td><td>55</td><td>2.7</td><td>2.0</td><td>0.0016</td></tr><tr><td>J0933-0819</td><td>-8.3</td><td>143.3</td><td>57411</td><td>533</td><td>3.1</td><td>2.0</td><td>0.0014</td><td>57128</td><td>697</td><td>2.9</td><td>2.0</td><td>0.0017</td></tr><tr><td>J0732+0150</td><td>1.8</td><td>113.1</td><td>55794</td><td>82</td><td>5.0</td><td>3.5</td><td>0.0010</td><td>55854</td><td>61</td><td>2.7</td><td>3.4</td><td>0.033</td></tr><tr><td>J0242+1101</td><td>11.0</td><td>40.6</td><td>56676</td><td>311</td><td>5.4</td><td>2.2</td><td>0.0060</td><td>56586</td><td>451</td><td>6.6</td><td>2.6</td><td>0.0021</td></tr><tr><td><b>J1826+1831</b></td><td><b>18.5</b></td><td><b>276.6</b></td><td>57672</td><td>151</td><td>2.9</td><td>2.5</td><td>0.0015</td><td><b>57636</b></td><td><b>178</b></td><td><b>3.0</b></td><td><b>2.5</b></td><td><b>0.0010</b></td></tr><tr><td>J1606+2717</td><td>27.3</td><td>241.7</td><td>58793</td><td>1</td><td>1.0</td><td>1.0</td><td>0.00076</td><td>58793</td><td>1</td><td>1.0</td><td>1.0</td><td>0.0017</td></tr><tr><td>J1800+3848</td><td>38.8</td><td>270.1</td><td>56590</td><td>3</td><td>1.7</td><td>2.4</td><td>0.0024</td><td>56590</td><td>3</td><td>1.9</td><td>2.6</td><td>0.0021</td></tr></tbody></table>

**Note.** List of radio-bright blazars for which a pretrial significance above $3\sigma$ 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 summarize the results of the search in terms of the best-fit central time of the flare $\hat{t}_0$, flare duration $\delta_t$, number of signal events $\hat{\mu}_s$, spectral index $\hat{\gamma}$, and pretrial $p$-value for the Gaussian-shaped and box-shaped signal time profiles. The most significant flare found assuming each of the considered time shapes is highlighted in bold.

with the highest weight. Indeed, as the algorithm searches for clustering of events in space and time, an ensemble of multiple events detected at close times rather than a single high-weight event is identified as a flare.

While none of the investigated sources alone is observed at a significant detection level when accounting for trials, a further study is performed to check for the presence of a cumulative excess. By performing background PEs in which the same source catalog is targeted, it is found, as shown in Figure 8, that the probability of finding 18 or more sources with a pretrial

significance greater than $3\sigma$ is 1.4% ($2.5\sigma$). This result provides an additional hint to the time-independent analysis, that a fraction of the blazars contained in the VLBI catalog could be neutrino emitters.

## 7. Multimessenger Flare 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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Sky map in equatorial coordinates showing the positions of the VLBI blazars coincident with the $3\sigma$ flares found in ANTARES data (orange dots). See Table 1 for a list of these flares. The two blazars with the most significant excesses are depicted in orange, while the red dot indicates the blazar studied in Section 7. The locations of the Galactic plane (dashed line) and Galactic center (star) are also shown.</figcaption></figure>

<table><caption>**Table 1**
Results of the Time-dependent Likelihood Scan</caption><thead><tr><th colspan="3">Source</th><th colspan="6">Results</th></tr><tr><th rowspan="2">Name</th><th colspan="2">δ</th><th colspan="5">Gaussian-shaped Time Profile</th><th colspan="5">Box-shaped Time Profile</th></tr><tr><th>(deg)</th><th>α</th><th>τ̂<sub>0</sub><br>(MJD)</th><th>δ̂<sub>t</sub><br>(days)</th><th>μ̂<sub>s</sub></th><th>γ̂</th><th>p-value</th><th>τ̂<sub>0</sub><br>(MJD)</th><th>δ̂<sub>t</sub><br>(days)</th><th>μ̂<sub>s</sub></th><th>γ̂</th><th>p-value</th></tr></thead><tbody><tr><td>J0112-6634</td><td>-66.6</td><td>18.1</td><td>58215</td><td>304</td><td>4.5</td><td>2.7</td><td>0.0026</td><td>58154</td><td>305</td><td>3.8</td><td>2.7</td><td>0.0097</td></tr><tr><td><b>J1355-6326</b></td><td><b>-63.4</b></td><td><b>208.9</b></td><td><b>56524</b></td><td><b>1041</b></td><td><b>7.9</b></td><td><b>2.8</b></td><td><b>0.00018</b></td><td>56091</td><td>905</td><td>6.0</td><td>2.9</td><td>0.0048</td></tr><tr><td>J0359-6154</td><td>-61.9</td><td>59.7</td><td>56316</td><td>78</td><td>5.4</td><td>3.5</td><td>0.0022</td><td>56321</td><td>112</td><td>5.7</td><td>3.5</td><td>0.0013</td></tr><tr><td>J0522-6107</td><td>-61.1</td><td>80.6</td><td>56221</td><td>42</td><td>4.6</td><td>3.4</td><td>0.0034</td><td>56232</td><td>59</td><td>4.9</td><td>3.4</td><td>0.0023</td></tr><tr><td>J1220-5604</td><td>-56.1</td><td>185.1</td><td>58406</td><td>18</td><td>2.8</td><td>2.6</td><td>0.00029</td><td>58413</td><td>22</td><td>0.4</td><td>2.2</td><td>0.0032</td></tr><tr><td>J1825-5230</td><td>-52.5</td><td>276.3</td><td>57265</td><td>600</td><td>5.4</td><td>2.7</td><td>0.0031</td><td>57188</td><td>959</td><td>5.5</td><td>2.7</td><td>0.0027</td></tr><tr><td>J0641-3554</td><td>-35.9</td><td>100.3</td><td>58084</td><td>16</td><td>2.9</td><td>3.0</td><td>0.0021</td><td>58081</td><td>19</td><td>3.0</td><td>3.0</td><td>0.0018</td></tr><tr><td>J1418-3509</td><td>-35.2</td><td>214.7</td><td>58120</td><td>11</td><td>3.3</td><td>2.9</td><td>0.0018</td><td>58121</td><td>14</td><td>2.9</td><td>2.8</td><td>0.0021</td></tr><tr><td>J1500-2358</td><td>-24.0</td><td>225.2</td><td>55846</td><td>4</td><td>3.7</td><td>2.3</td><td>0.0016</td><td>55847</td><td>6</td><td>3.7</td><td>2.2</td><td>0.0015</td></tr><tr><td>J0521-1737</td><td>-17.6</td><td>80.3</td><td>57332</td><td>1</td><td>2.0</td><td>1.9</td><td>0.0011</td><td>57333</td><td>1</td><td>2.0</td><td>1.9</td><td>0.0023</td></tr><tr><td>J2345-1555</td><td>-15.9</td><td>356.3</td><td>57653</td><td>460</td><td>3.2</td><td>2.6</td><td>0.0011</td><td>57784</td><td>404</td><td>2.4</td><td>2.7</td><td>0.0030</td></tr><tr><td>J1537-1259</td><td>-13.0</td><td>234.3</td><td>58201</td><td>46</td><td>2.6</td><td>2.0</td><td>0.0019</td><td>58201</td><td>55</td><td>2.7</td><td>2.0</td><td>0.0016</td></tr><tr><td>J0933-0819</td><td>-8.3</td><td>143.3</td><td>57411</td><td>533</td><td>3.1</td><td>2.0</td><td>0.0014</td><td>57128</td><td>697</td><td>2.9</td><td>2.0</td><td>0.0017</td></tr><tr><td>J0732+0150</td><td>1.8</td><td>113.1</td><td>55794</td><td>82</td><td>5.0</td><td>3.5</td><td>0.0010</td><td>55854</td><td>61</td><td>2.7</td><td>3.4</td><td>0.033</td></tr><tr><td>J0242+1101</td><td>11.0</td><td>40.6</td><td>56676</td><td>311</td><td>5.4</td><td>2.2</td><td>0.0060</td><td>56586</td><td>451</td><td>6.6</td><td>2.6</td><td>0.0021</td></tr><tr><td><b>J1826+1831</b></td><td><b>18.5</b></td><td><b>276.6</b></td><td>57672</td><td>151</td><td>2.9</td><td>2.5</td><td>0.0015</td><td><b>57636</b></td><td><b>178</b></td><td><b>3.0</b></td><td><b>2.5</b></td><td><b>0.0010</b></td></tr><tr><td>J1606+2717</td><td>27.3</td><td>241.7</td><td>58793</td><td>1</td><td>1.0</td><td>1.0</td><td>0.00076</td><td>58793</td><td>1</td><td>1.0</td><td>1.0</td><td>0.0017</td></tr><tr><td>J1800+3848</td><td>38.8</td><td>270.1</td><td>56590</td><td>3</td><td>1.7</td><td>2.4</td><td>0.0024</td><td>56590</td><td>3</td><td>1.9</td><td>2.6</td><td>0.0021</td></tr></tbody></table>

**Note.** List of radio-bright blazars for which a pretrial significance above $3\sigma$ 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 summarize the results of the search in terms of the best-fit central time of the flare $\hat{t}_0$, flare duration $\delta_t$, number of signal events $\hat{\mu}_s$, spectral index $\hat{\gamma}$, and pretrial $p$-value for the Gaussian-shaped and box-shaped signal time profiles. The most significant flare found assuming each of the considered time shapes is highlighted in bold.

with the highest weight. Indeed, as the algorithm searches for clustering of events in space and time, an ensemble of multiple events detected at close times rather than a single high-weight event is identified as a flare.

While none of the investigated sources alone is observed at a significant detection level when accounting for trials, a further study is performed to check for the presence of a cumulative excess. By performing background PEs in which the same source catalog is targeted, it is found, as shown in Figure 8, that the probability of finding 18 or more sources with a pretrial

significance greater than $3\sigma$ is 1.4% ($2.5\sigma$). This result provides an additional hint to the time-independent analysis, that a fraction of the blazars contained in the VLBI catalog could be neutrino emitters.

## 7. Multimessenger Flare 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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** 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 $\hat{\sigma}_t$ and $\hat{T}_0$ found in each case and reported in Table 1. Tracks (showers) are shown in blue (red).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Cumulative excess of $3\sigma$ 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 of observing 18 or more sources with a pretrial flare significance above $3\sigma$ is 1.4%.</figcaption></figure>

Observatory (OVRO; Richards et al. 2011) for those sources of Table 1 for which radio data are available. Only for the source J0242+1101 (PKS 0239+108) has a notable overlap in

time been found between the best-fit neutrino flare identified in this analysis 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 Large Area Telescope (LAT) γ-ray telescope and the IceCube neutrino telescope compatible with the direction of J0242+1101 have also been studied. Figure 9 also reports the adaptive binned γ-ray light curve, obtained from Fermi data using the method described in Lott et al. (2012) and Kramarenko et al. (2022). Remarkably, the most significant Fermi γ-ray flare for this blazar happened during the flaring emission detected in radio and the period 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 Pushkarev et al. (2010), Kramarenko et al. (2022). The time distribution of the IceCube tracklike events in the 10 yr point-source sample IceCube Collaboration (2021) with a 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 $\nu_{\mu}$-induced track with the notable high energy of 50 TeV was detected during the flare, its reported angular uncertainty radius being 1°4. Furthermore, it is worth reporting
<figure><img src="image_9.png" /><figcaption>**Figure 7.** 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 $\hat{\sigma}_t$ and $\hat{T}_0$ found in each case and reported in Table 1. Tracks (showers) are shown in blue (red).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Cumulative excess of $3\sigma$ 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 of observing 18 or more sources with a pretrial flare significance above $3\sigma$ is 1.4%.</figcaption></figure>

Observatory (OVRO; Richards et al. 2011) for those sources of Table 1 for which radio data are available. Only for the source J0242+1101 (PKS 0239+108) has a notable overlap in

time been found between the best-fit neutrino flare identified in this analysis 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 Large Area Telescope (LAT) γ-ray telescope and the IceCube neutrino telescope compatible with the direction of J0242+1101 have also been studied. Figure 9 also reports the adaptive binned γ-ray light curve, obtained from Fermi data using the method described in Lott et al. (2012) and Kramarenko et al. (2022). Remarkably, the most significant Fermi γ-ray flare for this blazar happened during the flaring emission detected in radio and the period 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 Pushkarev et al. (2010), Kramarenko et al. (2022). The time distribution of the IceCube tracklike events in the 10 yr point-source sample IceCube Collaboration (2021) with a 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\ \text{deg}^2$ are depicted. While there is no evidence of time clustering of the IceCube events, a $\nu_{\mu}$-induced track with the notable high energy of 50 TeV was detected during the flare, its reported angular uncertainty radius being $1^\circ 4$. Furthermore, it is worth reporting
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Multimessenger 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 tracklike (shower-like) events within 5° (10°) from this object. The box profile has been drawn using the best-fit values of $\delta_t$ and $\hat{T}_0$ found in this analysis. Second panel: OVRO radio light curve. Third panel: adaptive binned $\gamma$-ray light curve obtained from Fermi LAT data. Bottom panel: weighted time distribution of the IceCube tracklike 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.</figcaption></figure>

that three IceCube alerts Abbasi et al. (2023e), namely, 131165:9342044, 129933:32926212, and 128672:38561326, have a likelihood best-fit direction that 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-$\gamma$ flare coincidence in J0242+1101 is an a posteriori nonblind finding. It is still instructive to evaluate how likely such a correlation is to arise by chance. This study is conducted by means of PEs. In each PE, the time-dependent

search described in Section 6 is performed on a mock data set obtained by scrambling the real data in R.A. Each PE results in a list of fake best-fit neutrino flares, similar to Table 1. The chance probability of the correlation is then calculated by comparing such fake best-fit flares with the available OVRO and Fermi light curves as the percentage of flares fulfilling the following criteria:

1. the global maximum of their radio/$\gamma$ light curve falls within the ANTARES flare duration, $\hat{T}_0 \pm \delta_t$, and
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Multimessenger 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 tracklike (shower-like) events within 5° (10°) from this object. The box profile has been drawn using the best-fit values of $\delta_t$ and $\hat{T}_0$ found in this analysis. Second panel: OVRO radio light curve. Third panel: adaptive binned $\gamma$-ray light curve obtained from Fermi LAT data. Bottom panel: weighted time distribution of the IceCube tracklike 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.</figcaption></figure>

that three IceCube alerts Abbasi et al. (2023e), namely, 131165:9342044, 129933:32926212, and 128672:38561326, have a likelihood best-fit direction that 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-$\gamma$ flare coincidence in J0242+1101 is an a posteriori nonblind finding. It is still instructive to evaluate how likely such a correlation is to arise by chance. This study is conducted by means of PEs. In each PE, the time-dependent

search described in Section 6 is performed on a mock data set obtained by scrambling the real data in R.A. Each PE results in a list of fake best-fit neutrino flares, similar to Table 1. The chance probability of the correlation is then calculated by comparing such fake best-fit flares with the available OVRO and Fermi light curves as the percentage of flares fulfilling the following criteria:

1. the global maximum of their radio/$\gamma$ light curve falls within the ANTARES flare duration, $\hat{T}_0 \pm \delta_t$, and
2. 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 $\gamma$.

Thanks to the OVRO and Fermi observatories, both radio and $\gamma$-ray light curves are available for hundreds of blazars: the CGRaBS OVRO sample Hovatta et al. (2021) and weekly curves from the Fermi LAT light Curve Repository<fnref n="59" /> Abdollahi et al. (2023). 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 abovementioned characteristics; it appears in only $p=0.5\%$ of random realizations, 4% for either radio or $\gamma$-rays separately. This fraction of $p=0.5\%$ is the $p$-value, the chance coincidence probability of observing 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.

## 8. 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 tracklike neutrino candidates. The pair-counting method finds an excess of 59 neutrino–blazar pairs with a posttrial $p$-value of 0.03 ($2.2\sigma$). 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^{-\gamma}$, mildly constraining the blazar contribution to ~65% of the high-energy neutrino diffuse flux measured by the IceCube Collaboration and also supported by ANTARES data.

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\sigma$ (pretrial). The two most significant flares come from the blazars J1355–6326 (3.7$\sigma$ pretrial) and J1826+1831 (3.3$\sigma$ pretrial); see Table 1. However, neither is significant alone after correcting for multiple trials. A cumulative effect is detected with $p=0.014$ ($2.5\sigma$). 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 $\gamma$-ray light curves is $p=0.5\%$. This is a post hoc estimate that could be a hint of 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 blazars uniformly selected by their bright radio emission 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, and 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 be expected in the coming years. Radio observational programs are ongoing to aid in reliable identification of coincidences and better understanding of associations.

## Acknowledgments

We thank Eduardo Ros for helpful comments on the manuscript. This research has made use of data from the OVRO 40 m monitoring program, 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.

The authors acknowledge the financial support of the funding agencies: Centre National de la Recherche Scientifique (CNRS), Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Commission Européenne (FEDER fund and Marie Curie Program), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Région Alsace (contrat CPER), Région Provence-Alpes-Côte d'Azur, Département du Var and Ville de La Seyne-sur-Mer, France; Bundesministerium für Bildung und Forschung (BMBF), Germany; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; Executive Unit for Financing Higher Education, Research, Development and Innovation (UEFISCDI), Romania; grants PID2021-124591NB-C41, -C42, -C43 funded by MCIN/AEI/10.13039/501100011033 and, as appropriate, by "ERDF A way of making Europe," by the "European Union" or by the "European Union NextGenerationEU/PRTR," Programa de Planes Complementarios I+D+I (refs. ASFAE/2022/023, ASFAE/2022/014), Programa Prometeo (PROMETEO/2020/019) and GenT (refs. CIDEGENT/2018/034, /2019/043, /2020/049, /2021/23) of the Generalitat Valenciana, EU: MSC program (ref. 101025085), Spain; Ministry of Higher Education, Scientific Research and Innovation, Morocco; and the Arab Fund for Economic and Social Development, Kuwait. We also acknowledge the technical support of Ifremer, AIM, and Foselev Marine for the sea operation and the CC-IN2P3 for the computing facilities. Y.Y.K. acknowledges support from the M2FINDERS project, which has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (grant agreement No. 101018682). T.H. was supported by the Academy of Finland projects 317383, 320085, 322535, and 345899. S.K. acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 771282.

## ORCID iDs

J. Hofestädt https://orcid.org/0000-0002-7848-117X
G. Illuminati https://orcid.org/0000-0002-4138-8027
I. Kreykenbohm https://orcid.org/0000-0001-7335-1803
E. Leonora https://orcid.org/0000-0002-0536-3551

<footnote n="59">⁵⁹ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/</footnote>
2. 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 $\gamma$.

Thanks to the OVRO and Fermi observatories, both radio and $\gamma$-ray light curves are available for hundreds of blazars: the CGRaBS OVRO sample Hovatta et al. (2021) and weekly curves from the Fermi LAT light Curve Repository<fnref n="59" /> Abdollahi et al. (2023). 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 abovementioned characteristics; it appears in only $p=0.5\%$ of random realizations, 4% for either radio or $\gamma$-rays separately. This fraction of $p=0.5\%$ is the $p$-value, the chance coincidence probability of observing 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.

## 8. 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 tracklike neutrino candidates. The pair-counting method finds an excess of 59 neutrino–blazar pairs with a posttrial $p$-value of 0.03 ($2.2\sigma$). 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^{-\gamma}$, mildly constraining the blazar contribution to ~65% of the high-energy neutrino diffuse flux measured by the IceCube Collaboration and also supported by ANTARES data.

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\sigma$ (pretrial). The two most significant flares come from the blazars J1355–6326 (3.7$\sigma$ pretrial) and J1826+1831 (3.3$\sigma$ pretrial); see Table 1. However, neither is significant alone after correcting for multiple trials. A cumulative effect is detected with $p=0.014$ ($2.5\sigma$). 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 $\gamma$-ray light curves is $p=0.5\%$. This is a post hoc estimate that could be a hint of 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 blazars uniformly selected by their bright radio emission 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, and 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 be expected in the coming years. Radio observational programs are ongoing to aid in reliable identification of coincidences and better understanding of associations.

## Acknowledgments

We thank Eduardo Ros for helpful comments on the manuscript. This research has made use of data from the OVRO 40 m monitoring program, 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.

The authors acknowledge the financial support of the funding agencies: Centre National de la Recherche Scientifique (CNRS), Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Commission Européenne (FEDER fund and Marie Curie Program), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Région Alsace (contrat CPER), Région Provence-Alpes-Côte d'Azur, Département du Var and Ville de La Seyne-sur-Mer, France; Bundesministerium für Bildung und Forschung (BMBF), Germany; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; Executive Unit for Financing Higher Education, Research, Development and Innovation (UEFISCDI), Romania; grants PID2021-124591NB-C41, -C42, -C43 funded by MCIN/AEI/10.13039/501100011033 and, as appropriate, by “ERDF A way of making Europe,” by the “European Union” or by the “European Union NextGenerationEU/PRTR,” Programa de Planes Complementarios I+D+I (refs. ASFAE/2022/023, ASFAE/2022/014), Programa Prometeo (PROMETEO/2020/019) and GenT (refs. CIDEGENT/2018/034, /2019/043, /2020/049, /2021/23) of the Generalitat Valenciana, EU: MSC program (ref. 101025085), Spain; Ministry of Higher Education, Scientific Research and Innovation, Morocco; and the Arab Fund for Economic and Social Development, Kuwait. We also acknowledge the technical support of Ifremer, AIM, and Foselev Marine for the sea operation and the CC-IN2P3 for the computing facilities. Y.Y.K. acknowledges support from the M2FINDERS project, which has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (grant agreement No. 101018682). T.H. was supported by the Academy of Finland projects 317383, 320085, 322535, and 345899. S.K. acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 771282.

## ORCID iDs

J. Hofestädt https://orcid.org/0000-0002-7848-117X
G. Illuminati https://orcid.org/0000-0002-4138-8027
I. Kreykenbohm https://orcid.org/0000-0001-7335-1803
E. Leonora https://orcid.org/0000-0002-0536-3551

<footnote n="59">⁵⁹ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/</footnote>
P. Piattelli https://orcid.org/0000-0003-4748-6485

A. Plavin https://orcid.org/0000-0003-2914-8554

A. Pushkarev https://orcid.org/0000-0002-9702-2307

A. Sánchez-Losa https://orcid.org/0000-0001-9596-7078

M. Sanguineti https://orcid.org/0000-0002-7206-2097

F. Schüssler https://orcid.org/0000-0003-1500-6571

M. Spurio https://orcid.org/0000-0002-8698-3655

Th. Stolarczyk https://orcid.org/0000-0002-0551-7581

B. Vallage https://orcid.org/0000-0003-1255-8506

T. Hovatta https://orcid.org/0000-0002-2024-8199

S. Kiehlmann https://orcid.org/0000-0001-6314-9177

## References

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, PhRvL, 113, 101101

Aartsen, M. G., Ackermann, A., Adams, J., et al. 2018a, Sci, 361, 147

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, Sci, 361, eaat1378

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, ApJ, 886, 12

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020, PhRvL, 125, 121104

Abbasi, R., Ackermann, M., Adams, J., et al. 2021, PhRvD, 104, 022002

Abbasi, R., Ackermann, M., Adams, J., et al. 2022a, Sci, 378, 538

Abbasi, R., Ackermann, M., Adams, J., et al. 2022b, ApJ, 928, 50

Abbasi, R., Ackermann, M., Adams, J., et al. 2022c, ApJ, 928, 50

Abbasi, R., Ackermann, M., Adams, J., et al. 2023a, ApJ, 951, 45

Abbasi, R., Ackermann, M., Adams, J., et al. 2023b, Sci, 380, 1338

Abbasi, R., Ackermann, M., Adams, J., et al. 2023c, ApJS, 269, 25

Abbasi, R., Ackermann, M., Adams, J., et al. 2023d, ApJ, 954, 75

Abbasi, R., Ackermann, M., Adams, J., et al. 2023e, ApJ, 951, 45

Abdollahi, S., Ajello, M., Baldini, L., et al. 2023, ApJS, 265, 31

Adrian-Martinez, S., Al Samarai, I., Albert, A., et al. 2013, JCAP, 03, 006

Adrian-Martinez, S., Albert, A., André, M., et al. 2015, JCAP, 12, 014

Ageron, M., Aguilar, J.A., Al Samarai, I., et al. 2011, NIMPA, A656, 11

Aguilar, J. A., Albert, A., Ameli, F., et al. 2007, NIMPA, 570, 107

Aguilar, J. A., Al Samarai, I., Albert, A., et al. 2011, APh, 34, 652

Albert, A., André, M., Anghinolfi, M., et al. 2017a, PhRvD, 96, 062001

Albert, A., André, M., Anghinolfi, M., et al. 2017b, PhRvD, 96, 082001

Albert, A., André, M., Anghinolfi, M., et al. 2017c, AJ, 154, 275

Albert, A., André, M., Anghinolfi, M., et al. 2018a, ApJL, 868, L20

Albert, A., André, M., Anghinolfi, M., et al. 2018b, EPJC, 78, 669

Albert, A., André, M., Anghinolfi, M., et al. 2018c, ApJL, 853, L7

Albert, A., André, M., Anghinolfi, M., et al. 2019, ApJ, 879, 108

Albert, A., André, M., Anghinolfi, M., et al. 2021a, ApJ, 911, 48

Albert, A., André, M., Anghinolfi, M., et al. 2021b, JCAP, 01, 064

Albert, A., Alves, S., André, M., et al. 2023, PhLB, 841, 137951

Allakhverdyan, V. A., Avrorin, A. D., Avrorin, A. V., et al. 2023, PhRvD, 107, 042005

Amram, P., Anghinolfi, M., Anvar, S., et al. 2002, NIMPA, 484, 369

Beasley, A. J., Gordon, D., Peck, A. B., et al. 2002, ApJS, 141, 13

Berezinsky, V. 1977, in Proc. of the Int. Neutrino Conference, ed. H. Faissner, H. Reithler, & P. Zerwas (Braunschweig: Vieweg), 650

Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467

Boettcher, M. 2019, Galax. 7, 20

Boettcher, M. 2023, arXiv:2308.01083

Boettcher, M., Fu, M., Govenor, T., King, Q., & Roustazadeh, P. 2022, AcPPB, 15, 8

Braun, J., Baker, M., Dumm, J., et al. 2010, APh, 33, 175

Buson, S., Tramacere, A., Oswald, L., et al. 2023, arXiv:2305.11263

Buson, S., Tramacere, A., Pfeiffer, L., et al. 2022, ApJL, 933, L43

Eichler, D. 1979, ApJ, 232, 106

Fomalont, E. B., Petrov, L., MacMillan, D. S., et al. 2003, AJ, 126, 2562

Fusco, L. A., & Versari, F. 2019, ICRC (Madison, WI), 36, 891

Giommi, P., Glauch, T., Padovani, P., et al. 2020, MNRAS, 497, 865

Gordon, D., Jacobs, C., Beasley, A., et al. 2016, AJ, 151, 154

Hovatta, T., Lindfors, E., Kiehlmann, S., et al. 2021, A&A, 650, A83

IceCube Collaboration 2021, A&A, 650, A83

Kovalev, Y. Y., Petrov, L., Fomalont, E. B., & Gordon, D. 2007, AJ, 133, 1236

Kovalev, Y. Y., Plavin, A. V., & Troitsky, S. V. 2022, ApJL, 940, L41

Kramarenko, I. G., Pushkarev, A. B., Kovalev, Y. Y., et al. 2022, MNRAS, 510, 469

Kun, E., Bartos, I., Tjus Becker, J., et al. 2022, ApJ, 934, 180

Lott, B., Escande, L., Larsson, S., & Ballet, J. 2012, A&A, 544, A6

Murase, K. 2017, in Neutrino Astronomy: Current Status, Future Prospects, ed. T. Gaisser & A. Karle (Singapore: World Scientific), 15

Neronov, A., & Semikoz, D. V. 2016, APh, 75, 60

Petrov, L. 2011, AJ, 142, 105

Petrov, L. 2012, MNRAS, 419, 1097

Petrov, L. 2013, AJ, 146, 5

Petrov, L. 2021, AJ, 161, 14

Petrov, L., de Witt, A., Sadler, E. M., Phillips, C., & Horiuchi, S. 2019, MNRAS, 485, 88

Petrov, L., Gordon, D., Gipson, J., et al. 2009, JGeod, 83, 859

Petrov, L., Kovalev, Y. Y., Fomalont, E., & Gordon, D. 2005, AJ, 129, 1163

Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2006, AJ, 131, 1872

Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2008, AJ, 136, 580

Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2011a, MNRAS, 142, 35

Petrov, L., Phillips, C., Bertarini, A., et al. 2011b, MNRAS, 414, 2528

Piner, B. G., Pushkarev, A. B., Kovalev, Y. Y., et al. 2012, ApJ, 758, 84

Plavin, A., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2020, ApJ, 894, 101

Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2021, ApJ, 908, 157

Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2023, MNRAS, 523, 1799

Pushkarev, A. B., & Kovalev, Y. Y. 2012, A&A, 544, A34

Pushkarev, A. B., Kovalev, Y. Y., & Lister, M. L. 2010, ApJL, 722, L7

Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29

Schinzel, F. K., Petrov, L., Taylor, G. B., et al. 2015, ApJS, 217, 4

Shu, F., Petrov, L., Jiang, W., et al. 2017, ApJS, 230, 13

Zhou, B., Kamionkowski, M., & Liang, Y.-F. 2021, PhRvD, 103, 123018
P. Piattelli https://orcid.org/0000-0003-4748-6485

A. Plavin https://orcid.org/0000-0003-2914-8554

A. Pushkarev https://orcid.org/0000-0002-9702-2307

A. Sánchez-Losa https://orcid.org/0000-0001-9596-7078

M. Sanguineti https://orcid.org/0000-0002-7206-2097

F. Schüssler https://orcid.org/0000-0003-1500-6571

M. Spurio https://orcid.org/0000-0002-8698-3655

Th. Stolarczyk https://orcid.org/0000-0002-0551-7581

B. Vallage https://orcid.org/0000-0003-1255-8506

T. Hovatta https://orcid.org/0000-0002-2024-8199

S. Kiehlmann https://orcid.org/0000-0001-6314-9177

## References

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, PhRvL, 113, 101101

Aartsen, M. G., Ackermann, A., Adams, J., et al. 2018a, Sci, 361, 147

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, Sci, 361, eaat1378

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, ApJ, 886, 12

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020, PhRvL, 125, 121104

Abbasi, R., Ackermann, M., Adams, J., et al. 2021, PhRvD, 104, 022002

Abbasi, R., Ackermann, M., Adams, J., et al. 2022a, Sci, 378, 538

Abbasi, R., Ackermann, M., Adams, J., et al. 2022b, ApJ, 928, 50

Abbasi, R., Ackermann, M., Adams, J., et al. 2022c, ApJ, 928, 50

Abbasi, R., Ackermann, M., Adams, J., et al. 2023a, ApJ, 951, 45

Abbasi, R., Ackermann, M., Adams, J., et al. 2023b, Sci, 380, 1338

Abbasi, R., Ackermann, M., Adams, J., et al. 2023c, ApJS, 269, 25

Abbasi, R., Ackermann, M., Adams, J., et al. 2023d, ApJ, 954, 75

Abbasi, R., Ackermann, M., Adams, J., et al. 2023e, ApJ, 951, 45

Abdollahi, S., Ajello, M., Baldini, L., et al. 2023, ApJS, 265, 31

Adrian-Martinez, S., Al Samarai, I., Albert, A., et al. 2013, JCAP, 03, 006

Adrian-Martinez, S., Albert, A., André, M., et al. 2015, JCAP, 12, 014

Ageron, M., Aguilar, J.A., Al Samarai, I., et al. 2011, NIMPA, A656, 11

Aguilar, J. A., Albert, A., Ameli, F., et al. 2007, NIMPA, 570, 107

Aguilar, J. A., Al Samarai, I., Albert, A., et al. 2011, APh, 34, 652

Albert, A., André, M., Anghinolfi, M., et al. 2017a, PhRvD, 96, 062001

Albert, A., André, M., Anghinolfi, M., et al. 2017b, PhRvD, 96, 082001

Albert, A., André, M., Anghinolfi, M., et al. 2017c, AJ, 154, 275

Albert, A., André, M., Anghinolfi, M., et al. 2018a, ApJL, 868, L20

Albert, A., André, M., Anghinolfi, M., et al. 2018b, EPJC, 78, 669

Albert, A., André, M., Anghinolfi, M., et al. 2018c, ApJL, 853, L7

Albert, A., André, M., Anghinolfi, M., et al. 2019, ApJ, 879, 108

Albert, A., André, M., Anghinolfi, M., et al. 2021a, ApJ, 911, 48

Albert, A., André, M., Anghinolfi, M., et al. 2021b, JCAP, 01, 064

Albert, A., Alves, S., André, M., et al. 2023, PhLB, 841, 137951

Allakhverdyan, V. A., Avrorin, A. D., Avrorin, A. V., et al. 2023, PhRvD, 107, 042005

Amram, P., Anghinolfi, M., Anvar, S., et al. 2002, NIMPA, 484, 369

Beasley, A. J., Gordon, D., Peck, A. B., et al. 2002, ApJS, 141, 13

Berezinsky, V. 1977, in Proc. of the Int. Neutrino Conference, ed. H. Faissner, H. Reithler, & P. Zerwas (Braunschweig: Vieweg), 650

Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467

Boettcher, M. 2019, Galax. 7, 20

Boettcher, M. 2023, arXiv:2308.01083

Boettcher, M., Fu, M., Govenor, T., King, Q., & Roustazadeh, P. 2022, AcPPB, 15, 8

Braun, J., Baker, M., Dumm, J., et al. 2010, APh, 33, 175

Buson, S., Tramacere, A., Oswald, L., et al. 2023, arXiv:2305.11263

Buson, S., Tramacere, A., Pfeiffer, L., et al. 2022, ApJL, 933, L43

Eichler, D. 1979, ApJ, 232, 106

Fomalont, E. B., Petrov, L., MacMillan, D. S., et al. 2003, AJ, 126, 2562

Fusco, L. A., & Versari, F. 2019, ICRC (Madison, WI), 36, 891

Giommi, P., Glauch, T., Padovani, P., et al. 2020, MNRAS, 497, 865

Gordon, D., Jacobs, C., Beasley, A., et al. 2016, AJ, 151, 154

Hovatta, T., Lindfors, E., Kiehlmann, S., et al. 2021, A&A, 650, A83

IceCube Collaboration 2021, A&A, 650, A83

Kovalev, Y. Y., Petrov, L., Fomalont, E. B., & Gordon, D. 2007, AJ, 133, 1236

Kovalev, Y. Y., Plavin, A. V., & Troitsky, S. V. 2022, ApJL, 940, L41

Kramarenko, I. G., Pushkarev, A. B., Kovalev, Y. Y., et al. 2022, MNRAS, 510, 469

Kun, E., Bartos, I., Tjus Becker, J., et al. 2022, ApJ, 934, 180

Lott, B., Escande, L., Larsson, S., & Ballet, J. 2012, A&A, 544, A6

Murase, K. 2017, in Neutrino Astronomy: Current Status, Future Prospects, ed. T. Gaisser & A. Karle (Singapore: World Scientific), 15

Neronov, A., & Semikoz, D. V. 2016, APh, 75, 60

Petrov, L. 2011, AJ, 142, 105

Petrov, L. 2012, MNRAS, 419, 1097

Petrov, L. 2013, AJ, 146, 5

Petrov, L. 2021, AJ, 161, 14

Petrov, L., de Witt, A., Sadler, E. M., Phillips, C., & Horiuchi, S. 2019, MNRAS, 485, 88

Petrov, L., Gordon, D., Gipson, J., et al. 2009, JGeod, 83, 859

Petrov, L., Kovalev, Y. Y., Fomalont, E., & Gordon, D. 2005, AJ, 129, 1163

Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2006, AJ, 131, 1872

Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2008, AJ, 136, 580

Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2011a, MNRAS, 142, 35

Petrov, L., Phillips, C., Bertarini, A., et al. 2011b, MNRAS, 414, 2528

Piner, B. G., Pushkarev, A. B., Kovalev, Y. Y., et al. 2012, ApJ, 758, 84

Plavin, A., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2020, ApJ, 894, 101

Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2021, ApJ, 908, 157

Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2023, MNRAS, 523, 1799

Pushkarev, A. B., & Kovalev, Y. Y. 2012, A&A, 544, A34

Pushkarev, A. B., Kovalev, Y. Y., & Lister, M. L. 2010, ApJL, 722, L7

Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29

Schinzel, F. K., Petrov, L., Taylor, G. B., et al. 2015, ApJS, 217, 4

Shu, F., Petrov, L., Jiang, W., et al. 2017, ApJS, 230, 13

Zhou, B., Kamionkowski, M., & Liang, Y.-F. 2021, PhRvD, 103, 123018