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

CrossMark

# ANTARES Search for Point Sources of Neutrinos Using Astrophysical Catalogs: A Likelihood Analysis

A. Albert¹,², M. André³, M. Anghinolfi⁴, G. Anton⁵, M. Ardid⁶, J.-J. Aubert⁷, J. Aublin⁸, B. Baret⁸, S. Basa⁹, B. Belhorma¹⁰, V. Bertin⁷, S. Biagi¹¹, M. Bissinger⁵, J. Boumaaza¹², M. Bouta¹³, M. C. Bouwhuis¹⁴, H. Brânzaş¹⁵, R. Bruijn¹⁴,¹⁶, J. Brunner⁷, J. Busto⁷, A. Capone¹⁷,¹⁸, L. Caramete¹⁵, J. Carr⁷, V. Carretero¹⁹, S. Celli¹⁷,¹⁸, M. Chabab²⁰, T. N. Chau⁸, R. Cherkaoui El Moursli¹², T. Chiarusi²¹, M. Circella²², A. Coleiro⁸, M. Colomer-Molla⁸,¹⁹, R. Coniglione¹¹, P. Coyle⁷, A. Creusot⁸, A. F. Díaz²³, G. de Wasseige⁸, A. Deschamps²⁴, C. Distefano¹¹, I. Di Palma¹⁷,¹⁸, A. Domi⁴,²⁵, C. Donzaud⁸,²⁶, D. Dornic⁷, D. Drouhin¹,², T. Eberl⁵, N. El Khayati¹², A. Enzenhöfer⁷, P. Fermani¹⁷,¹⁸, G. Ferrara¹¹, F. Filippini²¹,²⁷, L. Fusco⁷,⁸, Y. Gatelet⁸, P. Gay⁸,²⁸, H. Glotin²⁹, R. Gozzini⁵,¹⁹, R. Gracía¹⁴, K. Graf⁵, C. Guidi⁴,²⁵, 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⁸,¹⁹, C. W. James³¹, B. Jisse-Jung¹⁴, M. de Jong¹⁴,³², P. de Jong¹⁴, M. Jongen¹⁴, M. Kadler³³, O. Kalekin⁵, U. Katz⁵, N. R. Khan-Chowdhury¹⁹, A. Kouchner⁸, I. Kreykenbohm³⁴, V. Kulikovskiy⁴,³⁵, R. Lahmann⁵, R. Le Breton⁸, D. Lefèvre³⁶, E. Leonora³⁷, G. Levi²¹,²⁷, M. Lincetto⁷, D. Lopez-Coto³⁸, S. Loucatos⁸,³⁹, L. Maderer⁸, J. Manczak¹⁹, M. Marcelin⁹, A. Margiotta²¹,²⁷, A. Marinelli⁴⁰, J. A. Martínez-Mora⁶, S. Mazzou²⁰, K. Melis¹⁴,¹⁶, P. Migliozzi⁴⁰, M. Moser⁵, A. Moussa¹³, R. Muller¹⁴, L. Nauta¹⁴, S. Navas³⁸, E. Nezri⁹, A. Nuñez-Castiñeyra⁷,⁹, B. O’Fearraigh¹⁴, M. Organokov¹, G. E. Păvălas¹⁴,¹⁵, C. Pellegrino²¹,⁴¹,⁴², M. Perrin-Terrin⁷, P. Piattelli¹¹, C. Pieterse¹⁹, C. Poirè⁶, V. Popa¹⁵, T. Pradier¹, N. Randazzo³⁷, S. Reck⁵, G. Riccobene¹¹, F. Salesa Greus¹⁹, D. F. E. Samtleben¹⁴,³², A. Sánchez-Losa²², M. Sanguineti⁴,²⁵, P. Sapienza¹¹, J. Schnabel⁵, F. Schüssler³⁹, M. Spurio²¹,²⁷, Th. Stolarczyk³⁹, M. Taiuti⁴,²⁵, Y. Tayalati¹², T. Thakore¹⁹, S. J. Tingay³¹, B. Vallage⁸,³⁹, V. Van Elewyck⁸,⁴³, F. Versari⁸,²¹,²⁷, S. Viola¹¹, D. Vivolo⁴⁰,⁴⁴, J. Wilms³⁴, A. Zegarelli¹⁷,¹⁸, J. D. Zornoza¹⁹, J. Zúñiga¹⁹

(ANTARES Collaboration), and

S. Buson⁴⁵

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

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

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

⁴ INFN—Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

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

⁶ 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é de Paris, CNRS, Astroparticule et Cosmologie, F-75006 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

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

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

¹³ 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 Roma, P. le Aldo Moro 2, 00185 Roma, Italy

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

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

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

²¹ INFN—Sezione di Bologna, Viale Berti-Pichat 6/2, 40127 Bologna, Italy

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

²³ Department of Computer Architecture and Technology/CITIC, University of Granada, 18071 Granada, Spain

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

²⁵ Dipartimento di Fisica dell'Università, Via Dodecaneso 33, 16146 Genova, Italy

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

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

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

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

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

³¹ International Centre for Radio Astronomy Research—Curtin University, Bentley, WA 6102, Australia

³² 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

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

³⁵ Moscow State University, Skobeltsyn Institute of Nuclear Physics, Leninskie gory, 119991 Moscow, Russia

³⁶ 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

³⁷ 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, Via Cintia 80126 Napoli, Italy

⁴¹ Museo Storico della Fisica e Centro Studi e Ricerche Enrico Fermi, Piazza del Viminale 1, 00184, Roma, Italy

⁴² INFN—CNAF, Viale C. Berti Pichat 6/2, 40127, Bologna, Italy
Abstract

A search for astrophysical pointlike neutrino sources using the data collected by the ANTARES detector between 2007 January 29 and 2017 December 31 is presented. A likelihood method is used to assess the significance of an excess of muon neutrinos inducing track-like events in correlation with the location of a list of possible sources. Different sets of objects are tested in the analysis: (a) a subsample of the Fermi 3LAC catalog of blazars, (b) a jet-obscured population of active galactic nuclei, (c) a sample of hard X-ray selected radio galaxies, (d) a star-forming galaxy catalog, and (e) a public sample of 56 very-high-energy track events from the IceCube experiment. None of the tested sources shows a significant association with the sample of neutrinos detected by ANTARES. The smallest *p*-value is obtained for the catalog of radio galaxies with an equal-weights hypothesis, with a pre-trial *p*-value equivalent to a 2.8σ excess, which is equivalent to 1.6σ post-trial. In addition, the results of a dedicated analysis for the blazar MGJ3J225517+2409 are also reported: this source is found to be the most significant within the Fermi 3LAC sample, with five ANTARES events located less than one degree from the source. This blazar showed evidence of flaring activity in Fermi data, in spacetime coincidence with a high-energy track detected by IceCube. An a posteriori significance of 2.6σ for the combination of ANTARES and IceCube data is reported.

*Unified Astronomy Thesaurus concepts: Neutrino astronomy (1100); Particle astrophysics (96); Astronomy data analysis (1858); Catalogs (205)*

1. Introduction

The origin and nature of cosmic rays at very high energy is a long-standing puzzle that is still not completely resolved after decades of theoretical and experimental efforts. The fact that cosmic rays are mostly charged particles, experiencing magnetic deflections during their propagation, makes the identification of their sources very difficult.

In contrast, neutrinos (ν) are neutral particles that can travel without absorption or deflection from their source to reach the Earth. Neutrinos of cosmic origin are expected to be produced via the decay of charged pions and kaons, generated in hadronic interactions of cosmic rays with gas or radiation at their acceleration sites or during their propagation. The decay of neutral pions also produces gamma rays (γ) that should be detectable, thus intimately connecting these three messengers: cosmic rays, neutrinos, and gamma rays.

The observation of a diffuse astrophysical neutrino flux was established in 2013 by the IceCube Collaboration with a high level of significance (Aartsen et al. 2014), with the detection of an excess of neutrinos above the expected background. Their spatial distribution compatible with isotropy favors an extra-galactic origin.

Moreover, the first association between a high-energy neutrino and a cosmic source was reported by IceCube (Aartsen et al. 2018a) in 2017 September, when a through-going muon track with a deposited energy of ~24 TeV was detected in spatial coincidence with the position of the blazar TXS 0506+056. The high-energy neutrino occurred during a period of intense gamma-ray emission observed by the Fermi-LAT and MAGIC telescopes (Aartsen et al. 2018a). A subsequent study (Aartsen et al. 2018b) based on IceCube's archival data in the direction of the blazar reported the presence of a 3.5σ candidate TeV neutrino flare between 2014 September and 2015 March positionally consistent with TXS 0506+056. ANTARES data do not show strong evidence of neutrinos in the direction of this blazar (Albert et al. 2018).

However, during the period of the IceCube neutrino flare, the gamma-ray flux of the source was one order of magnitude lower than its value during the period of intense gamma-ray

activity and it showed no sign of time variability. If the neutrino signal observed in the direction of TXS 0506+056 is truly of astrophysical origin and related to the blazar, this implies that the relation between the observed gamma-ray flux and the high-energy neutrino flux is not simple, and can depend on the individual properties of each source. More multi-messenger observations are then needed to constrain the models and gather evidence to identify the sources of high-energy neutrinos.

The recent findings using IceCube data (Padovani et al. 2016; Garrappa & Buson 2019; Giommi et al. 2020) support the idea that blazars could produce an observable high-energy neutrino flux. Among the various potential sources of high-energy neutrinos, active galactic nuclei (AGNs), and more specifically blazars, are therefore of particular interest. Indeed, when one of the relativistic jets emitted by the central black hole is pointing close to Earth's line of sight, the relativistic Doppler effect can produce a boosted flux of electromagnetic radiation and potentially related neutrinos if hadronic interactions occur within the jet. The mechanism of production of high-energy neutrinos from blazars has been studied in numerous papers (Mannheim 1995; Halzen & Zas 1997; Mücke et al. 2003; Murase et al. 2014; Tavecchio & Ghisellini 2015) and is still an active field of investigation.

From the experimental point of view, different methods have been used to search for neutrino point sources on the full observable sky (Albert et al. 2017; Albert & Aartsen et al. 2020; Aartsen et al. 2020a), including searches over a predefined target list (Aartsen et al. 2017; Huber & Krings 2017; Huber 2019). The large number of trials associated with the full-sky searches implies that only very bright objects can be detected with >5σ significance with such methods. Searching among a limited set of directions, for example from a catalog of sources, is an efficient way to reduce the trial factors and also allows population studies to be performed.

The present analysis uses a likelihood stacking method, allowing the significance of an excess of neutrino events in the direction of a selected sample of sources to be assessed. The potential neutrino emitters are considered together, as a

<footnote n="43"><sup>43</sup> Institut Universitaire de France, 75005 Paris, France</footnote>

<footnote n="44"><sup>44</sup> Dipartimento di Fisica dell'Università Federico II di Napoli, Via Cintia 80126, Napoli, Italy</footnote>

<footnote n="45"><sup>45</sup> Universität Würzburg, Institut für Theoretische Physik und Astrophysik, Campus Hubland Nord, Emil-Fischer-Str. 31 97074 Würzburg, Germany

*Received 2020 December 29; revised 2021 February 3; accepted 2021 February 9; published 2021 April 14*</footnote>
common population from which neutrinos could be produced. The ANTARES data provide an independent measurement with respect to the IceCube data, offering a complementary sensitivity in the Southern Hemisphere (see Figure 1 in Albert & Aartsen et al. 2020).

The ANTARES data collected between 2007 and 2017 are used to find possible correlations with several catalogs. The Fermi 3LAC catalog (Ackermann et al. 2015) is used to test the blazar origin hypothesis; other catalogs of astrophysical objects are also considered: dust-obscured AGNs (Maggi et al. 2016), very bright radio galaxies (Bassani et al. 2016), and star-forming galaxies (Ackermann et al. 2012). As an additional and independent test, the correlation of ANTARES neutrinos with 51 track events from the IceCube public high-energy neutrino candidates (Haack & Wiebusch 2017; Kopper 2017) is also evaluated.

The paper is organized as follows: the data set and the method are presented in Sections 2 and 3, the target catalogs are described in Section 4, and the results are shown in Section 5. An investigation of the most interesting individual objects is presented in subsection 5.1, and a dedicated analysis of the blazar MG3 J225517+2409 is performed in subsection 5.2. Finally, the results are summarized and discussed in Section 6.

## 2. Data Set

The ANTARES neutrino telescope is a water Cherenkov detector operating since 2007 on the bottom of the Mediterranean Sea, 40 km offshore from Toulon (France), which detects high-energy neutrinos by observing the passage in water of relativistic particles produced by interactions of neutrinos in the vicinity of the instrumented volume (Ageron et al. 2011).

The data sample used in the present analysis is presented in Albert et al. (2017) and consists of 8754 track-like events detected by ANTARES between 2007 January 29 and 2017 December 31, for an overall livetime of 3125.4 days. These events are almost 100% induced by charged-current interactions of $\nu_\mu$, the detected track corresponding to the path of the emitted muon.

The search for pointlike sources requires a very good angular resolution. The event selection has been optimized to detect a pointlike source with an energy spectrum $\propto E^{-2}$. The present data selection uses the same cuts that were defined in the 9 yr search (Albert et al. 2017), applied to an extended data set that includes two additional years. An in-depth description of the selection criteria and data/Monte Carlo comparison can be found in Albert et al. (2017), so only the most relevant parameters are described in the following.

The candidate neutrinos are selected by first considering mainly up-going directions, where the reconstructed zenith angle $\theta$ is reduced to values $\cos\theta > -0.1$, and by further cutting on the estimated angular error $\beta$ and the quality parameter $\Lambda$ of the likelihood reconstruction. The atmospheric muons, which constitute the dominant contribution to the background, are mostly suppressed by the cut on the quality parameter $\Lambda$ (~4 orders of magnitude) and by that on zenith angle (~1 additional order of magnitude).

After the selection procedure, a residual contamination of ~13% is estimated (Albert et al. 2017) to come from mis-reconstructed atmospheric muons. In comparison, about ~26 signal events would be expected in total from a reference

<figure><img src="image_2.png" /><figcaption>Figure 1. Cumulative probability to reconstruct the direction of a track event within an angle $\alpha$ with respect to the true Monte Carlo neutrino direction. The thick solid black line corresponds to the average value for an $E^{-2}$ energy spectrum, while the thinner lines are obtained for different energy ranges (still with $E^{-2}$ weights included).</figcaption></figure>

diffuse astrophysical $E^{-2}$ flux of $\nu_\mu + \bar{\nu}_\mu$ with a normalization of $10^{-8} \text{ GeV}^{-1} \text{ cm}^{-2} \text{ s}^{-1}$ at 1 GeV.

The energy estimator used in the analysis is based on the measure of the energy deposited by the muon in water (Schnabel 2013). The selected tracks have estimated neutrino energies between ~100 GeV and ~1 PeV, with a median angular resolution ranging from 1° at energies below 1 TeV to better than 0.4° above 10 TeV (see Figure 1).

## 3. Description of the Method

### 3.1. Test Statistic

An extended maximum likelihood method (Barlow 1990) is used to reject the null hypothesis $H_0$ where only background is assumed to be present, from an alternative hypothesis $H_1$ of signal plus background. The log-likelihoods for both hypotheses $H_0$ and $H_1$ are written respectively as

$$ 
\begin{aligned}
\ln \mathcal{L}(x|H_0) &= \sum_{i}^{N} \ln [\mu_b B(x_i)] - \mu_b \\
\ln \mathcal{L}(x|H_1) &= \sum_{i}^{N} \ln [\mu_s S(x_i) + \mu_b B(x_i)] - \mu_s - \mu_b,
\end{aligned}
\quad (1) $$

where $N$ is the total number of observed neutrino candidates, $S(x_i)$ is the probability density function (PDF) of the signal, and $B(x_i)$ the PDF of the background. The free parameters are the estimated numbers of signal $\mu_s$ and background $\mu_b$ events. The test statistic is then defined as

$$ \lambda = \ln \left( \frac{\max(\mathcal{L}(x|H_1))}{\max(\mathcal{L}(x|H_0))} \right). \quad (2) $$

### 3.2. Probability Density Functions of Signal and Background

#### 3.2.1. Signal

For a given catalog, the signal term is proportional to the expected number of neutrinos that should be detected as track-like events by ANTARES, from the directions of the $N_{\text{sources}}$ objects, each one of them emitting a $\nu_\mu + \bar{\nu}_\mu$ neutrino flux
$\Phi_j(E) = \Phi_j^0(\frac{E}{1 \text{ GeV}})^{-\gamma}$, the flux normalization $\Phi_j^0$ (GeV$^{-1}$ cm$^{-2}$ s$^{-1}$) being different for each source j. The spectral index $\gamma$ is assumed to be $\gamma = 2.0$ in the analysis, and is fixed in the likelihood maximization.

The signal PDF is then written as

$$S(x_i) = \frac{1}{\sum w_j} \sum_{j=1}^{N_{\text{sources}}} w_j s_j(x_i), \qquad (3)$$

where $s_j(x_i)$ is the contribution of the jth object, with weight $w_j$. The weight of the jth source is defined as

$$w_j = w_j^{\text{model}} \mathcal{A}(\delta_j), \qquad (4)$$

where $\mathcal{A}(\delta)$ is the declination-dependent acceptance of the neutrino track sample computed for an $E^{-2}$ energy spectrum.

The function $\mathcal{A}(\delta)$ allows the computation of the number of signal events expected for a single point source having a flux normalization $\Phi_0$ simply via the relation

$$n(\delta) = \Phi_0 \mathcal{A}(\delta). \qquad (5)$$

For a given energy spectrum $\Phi(E) = \Phi_0(\frac{E}{1 \text{ GeV}})^{-\gamma}$, the acceptance is computed from the effective area (Albert et al. 2017) of the detector, $\mathcal{A}_{\text{eff}}(E, \delta, t)$, with the relation

$$\mathcal{A}(\delta) = \iint dt \, dE \, \mathcal{A}_{\text{eff}}(E, \delta, t) \left(\frac{E}{1 \text{ GeV}}\right)^{-\gamma}, \qquad (6)$$

where the integral is performed over a sufficiently large energy range (100 GeV–10 PeV), and for the total lifetime considered in the analysis (3125.4 days). The acceptance obtained for a pure $E^{-2}$ energy spectrum is shown in Figure 5 of Albert et al. (2017).

The weighting scheme aims to set the weights to be proportional to the true neutrino flux. Since this is unknown, two different weighting assumptions are considered:

1. $w_j^{\text{model}} = \Phi_j^0$

2. $w_j^{\text{model}} = 1$ (equal-flux assumption),

where $\Phi_j^0$ is the photon flux observed in a relevant part of the electromagnetic spectrum, to be used as a proxy for the neutrino flux. As the relation between the flux emitted in a given band of the electromagnetic spectrum and the neutrino flux is rather uncertain, the equal-flux assumption ($w_j^{\text{model}} = 1$) provides an additional model-independent test, which could be more sensitive if the weights in the likelihood are very different from the true (unknown) ones.

The individual source contribution to the signal for the ith event is

$$s_j(x_i) = \mathcal{F}(\alpha_{ij}, E_i, \beta_i), \qquad (7)$$

where $\mathcal{F}$ is the energy-dependent point-spread function, defined as the probability density for the reconstructed event direction to lie within an angular distance $\alpha_{ij}$ from the true source direction.

Figure 1 shows that on average, for an $E^{-2}$ spectrum, 50% of the track events in the point-source sample are reconstructed within an angle of 0.4° of the true neutrino direction.

The most relevant parameters that determine the shape of the point-spread function are the estimated energy $E_i$ and the angular reconstruction uncertainty $\beta_i$ coming from the track reconstruction procedure.

The function $\mathcal{F}$ is evaluated using a collection of 3D histograms built from Monte Carlo simulations, in declination

bands of $\Delta \sin \delta = 0.2$. For a given track event with declination $\delta$, energy $E$, and angular uncertainty $\beta$, the histogram corresponding to the correct declination band is used, and the value of the signal is obtained by linear interpolation in the $(\alpha, E, \beta)$ space.

### 3.2.2. Background

A conservative approach is used to determine the background term. As the contribution of an astrophysical signal is expected to be small, real data are used to build the background PDF. The limited statistics available in real data as opposed to Monte Carlo do not allow the construction of a three-dimensional histogram, therefore the expression of the PDF is factorized as

$$B(x_i) = \mathcal{B}(\sin \delta_i) f_b(E_i, \beta_i), \qquad (8)$$

where $\mathcal{B}(\sin \delta)$ is the observed distribution of the sine of the declination, and $f_b$ the distribution of energy $E$ and angular uncertainty $\beta$ derived in data. The distribution of ANTARES track events measured during more than ten years is independent of the right ascension (see Figure 2 top left), because of the Earth’s rotation averaging the nonuniformity of the detector’s exposure in local coordinates.

Figure 2 (bottom left) shows the rate of events as a function of $\sin \delta$, fitted by a polynomial function that is used in the likelihood computation. The value of the function $f_b(E, \beta)$ is evaluated by linear interpolation in the 2D histogram shown in Figure 2 (right).

## 4. Target Sources

The different catalogs used in the analysis are presented in the following. Figure 3 shows the position of all the considered astrophysical objects on a sky map in equatorial coordinates, together with the positions of the candidate neutrinos detected by ANTARES.

### 4.1. Fermi 3LAC Blazars

The Fermi 3LAC catalog (Ackermann et al. 2015) contains blazars detected in gamma rays in the range 1–300 GeV by the Fermi-LAT in four years of operation. The previous release (Fermi 2LAC) has been used by the IceCube Collaboration in a stacking analysis (Aartsen et al. 2020a), where several subclasses of objects have been tested. A similar approach is followed here: in addition to the full Fermi 3LAC catalog, two additional subsamples are also considered, as described below.

Blazars have a typical spectral energy distribution (SED) of their electromagnetic spectrum that presents two bumps: a low-energy peak between IR and UV–X due to synchrotron emission of electrons, and a high-energy peak in the hard X-ray/gamma-ray band that can be produced by both hadronic and leptonic processes. Blazars are classified into flat-spectrum radio quasars (FSRQs) and BL Lacs, according to the presence (FSRQs) or absence (BL Lacs) of optical broad emission lines on top of the continuum. When this classification is not possible, a blazar is labelled as a blazar candidate of uncertain type (BCU) in the catalog.

In addition, a second classification can be made using the position of the synchrotron peak to distinguish low synchrotron peak (LSP) from intermediate and high synchrotron peaks (ISP and HSP respectively). This classification was proposed as an extension of the LBL/HBL classification of BL Lacs originally introduced by Padovani & Giommi (1995).
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Distribution of the right ascension (RA) and of $\sin\delta$ (left plots) obtained for events passing the analysis cuts in the data collected by ANTARES between 2007 January and 2017 December. The solid line represents the polynomial fit function used in the likelihood evaluation. The two-dimensional histogram representing the function $f_b(E, \beta)$ is shown in color code on the right plot (arbitrary units).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Sky map in equatorial coordinates of the 8754 ANTARES track events (small black circles), together with the considered target sources: Fermi 3LAC blazars (blue circles), radio galaxies (orange stars), star-forming galaxies (violet triangles), dust-obscured AGNs (red squares), and IceCube tracks (red empty circles). For astrophysical sources, the radius of each marker is proportional to the square root of the weight used in the likelihood. For IceCube tracks, the radius of each circle is proportional to the angular error (displayed angular size not to scale on the figure for clarity). The Galactic center is indicated by a black star and the Galactic plane by a blue dashed line. The blazar MG3 J225517+2409 and the radio galaxy 3C 403 have their positions indicated by their names. The blue shaded area corresponds to declinations out of the field of view of ANTARES.</figcaption></figure>

In addition, the objects that are flagged as doubtful in the Fermi catalog are removed from the analysis (see Ackermann et al. 2015, Section 2). The different (sub)samples that could be used in the analysis are:

1. Fermi 3LAC: 1420 objects, with 1255 in the field of view (FoV) of ANTARES,

2. FSRQ: 414 objects (382 in FoV),

3. BL Lacs: 604 objects (509 in FoV),

4. LSP: 666 objects (604 in FoV),

5. ISP+HSP: 685 objects (590 in FoV).

The overlap between subsamples is indicated in Table 1.

As can be seen, most of the FSRQs are classified as LSP and the majority of BL Lacs are classified as ISP+HSP. In order to reduce the number of trials, the analysis is then restricted to the

FSRQ and BL Lac subsamples only, in addition to the full Fermi 3LAC sample.

For the implementation in the likelihood, the measured slope $\gamma$ of the gamma-ray energy spectrum is taken into account in the weights. This is done by using the energy flux: $w^{\text{model}} = \int E \times \Phi_0(E/E_0)^{-\gamma} dE$ where the integral runs over the energy range 1–100 GeV.

## 4.2. Dust-obscured AGNs

The authors of Maggi et al. (2016) suggest considering a specific class of AGNs where the jet axis points toward the Earth but passes through a large quantity of dust or gas, thus potentially enhancing the neutrino production rate via hadronic p–p interactions.
<table><caption>**Table 1**
Repartition of the Fermi 3LAC Objects into the Different Subsamples</caption><thead><tr><th></th><th>FSRQ</th><th>BL Lac</th><th>BCU</th><th>Total</th></tr></thead><tbody><tr><td>LSP</td><td>365</td><td>148</td><td>153</td><td>666</td></tr><tr><td>ISP+HSP</td><td>44</td><td>440</td><td>201</td><td>685</td></tr><tr><td>Unclassified</td><td>5</td><td>16</td><td>48</td><td>69</td></tr><tr><td>Total</td><td>414</td><td>604</td><td>402</td><td>1420</td></tr></tbody></table>

The sample is built by cross-correlating objects from a catalog of radio-emitting galaxies (van Velzen et al. 2012) with the Fermi 2LAC catalog (Ackermann et al. 2011). The morphologies that are incompatible with a jet viewed face-on are rejected. The final selection criterion uses the observed X-ray flux of the objects, which should be low compared to the radio flux because of the attenuation in the surrounding matter.

The 25% weakest X-ray emitters are then selected, leading to a small catalog of 15 objects, containing three FSRQs, one ultraluminous IR galaxy, and 11 BL Lacs.

To take into account the distance and the intrinsic luminosity, the radio flux density measured at 1.4 GHz is used as a weight in the likelihood: $w^{\text{model}} = \Phi_{\text{Radio}}$.

### 4.3. Radio Galaxies

Radio galaxies are a subclass of active galaxies that could be considered as possible sources of high-energy neutrinos (Hooper 2016). The neutrinos could be produced by the hadronic interaction of accelerated cosmic rays within the jets or in the giant lobes at the ends of the jets.

Radio galaxies have an active nucleus emitting jets of relativistic plasma, but contrary to blazars, the direction of the jet axis is not close to our line of sight. Radio galaxies are then in general less luminous than blazars individually, but they are more abundant in the local universe.

A sample of 65 hard X-ray selected radio galaxies suggested by L. Bassani is considered in this study. The authors of Bassani et al. (2016) looked at extended radio galaxies in the sample of soft gamma-ray (20–100 keV) selected AGNs in either the INTEGRAL or Swift-BAT survey. The selection procedure yields a sample that contains the brightest and most accretion-efficient radio galaxies in the local sky, with very powerful jets producing an extended double-lobed radio morphology.

Among the published list of objects, one source without redshift information is removed, and Centaurus A is excluded because it cannot be considered as a point source (the angular extent of its radio lobes is ~10°).

The weight $w^{\text{model}} = L_{\gamma}/d^2$ is used, where $L_{\gamma}$ is the luminosity measured in soft gamma rays by Swift-BAT (except for five galaxies where the luminosity measured by INTEGRAL is used) and $d$ is the luminosity distance to the source. The luminosity distance is computed for each galaxy adopting the redshift indicated in Bassani et al. (2016), using a standard flat $\Lambda$CDM cosmology with parameters $H_0 = 71$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.27$, and $\Omega_{\Lambda} = 0.73$ (Hinshaw et al. 2013).

### 4.4. Star-forming Galaxies

A catalog of 64 star-forming galaxies (SFGs) published by the Fermi-LAT Collaboration (Ackermann et al. 2012) was

used recently by the Pierre Auger Collaboration, reporting a 3.9σ correlation with ultrahigh-energy cosmic rays with energy above 39 EeV (Aab et al. 2018). Each SFG does not necessarily have an active nucleus, but has a strong star formation rate traced by the intense infrared radiation. The dense interstellar gas can act as target material for the cosmic rays that are trapped in the galactic magnetic field, thus potentially producing neutrinos.

The dominant galaxies in terms of emitted infrared flux are M82, NGC 253, and NGC 4945. Due to the selection procedure (detection in gamma rays by Fermi-LAT), most of the objects present in the catalog are very nearby: five objects are located within a radius of 5 Mpc, 14 galaxies are within 10 Mpc, and 51 within 100 Mpc. A caveat of the use of this catalog is then its incompleteness due to the gamma-ray selection.

The following weights are used in the likelihood: $w^{\text{model}} = L_{\text{IR}}/d^2$ where $L_{\text{IR}}$ is the total IR luminosity (8–1000 µm) and $d$ is the luminosity distance to the source.

### 4.5. IceCube High-energy Tracks

The high-energy tracks measured by IceCube are an interesting sample for a stacking analysis: they are likely to be of astrophysical origin and have for the majority a median angular resolution of ~1°. Their directions are then in general a sufficiently good proxy for the location of neutrino sources.

The sample consists of 56 events, 55 of which are in the FoV of ANTARES: 35 tracks from the IceCube 8 yr sample of up-going muons with deposited energy above 200 TeV (Haack & Wiebusch 2017) and 21 additional tracks from the 6 yr HESE sample (Kopper 2017). The majority of those events are located in the Northern Hemisphere, and have been reconstructed with angular errors varying from ~0.5° up to ~5°. A complementary study of the correlation between IceCube high-energy tracks and ANTARES point-source data is reported in Illuminati et al. (2019), where each track is considered separately.

The IceCube candidate neutrinos have individual angular errors that need to be accounted for when searching for point-source excesses. The likelihood method is then modified as follows: an independent fit is performed for each IceCube track in the sample, where the location in (R.A., δ) is let free to vary in the fit within a range of ±2σ separately for each coordinate, where σ is the angular error reported by IceCube (for the 8 yr up-going muon sample, the 50% containment value reported in the article is multiplied by a factor 1.5).

The test statistic is then defined as the sum over all the individual test statistics:

$$ \lambda = \frac{1}{\sum_j w_j} \sum_{j=1}^{N_{\text{sources}}} w_j \lambda_j, \qquad (9) $$

where $\lambda_j$ is the local test statistic obtained for the jth IceCube track, carrying a weight $w_j = A(\delta_j)$ equal to the ANTARES acceptance at the corresponding declination.

The global significance is evaluated via the total test statistic $\lambda$; nevertheless, the relative contribution of each IceCube track can be examined by looking at the distribution of individual test statistics $\lambda_j$.
<table><caption>**Table 2**
Results of the Likelihood Stacking Analysis</caption><thead><tr><th rowspan="2">Catalog</th><th colspan="4">Equal Weighting</th><th colspan="4">Flux Weighting</th></tr><tr><th>λ</th><th>p</th><th>P</th><th>Φ<sup>UL</sup><sub>90%</sub></th><th>λ</th><th>p</th><th>P</th><th>Φ<sup>UL</sup><sub>90%</sub></th></tr></thead><tbody><tr><td>Fermi 3LAC All Blazars</td><td>6.1</td><td>0.19</td><td>0.83</td><td>4.3</td><td>0.21</td><td>0.85</td><td>1.0</td><td>2.1</td></tr><tr><td>Fermi 3LAC FSRQs</td><td>0.83</td><td>0.57</td><td>0.97</td><td>2.2</td><td>~0</td><td>~1</td><td>1.0</td><td>1.8</td></tr><tr><td>Fermi 3LAC BL Lacs</td><td>8.3</td><td>0.088</td><td>0.64</td><td>4.8</td><td>0.84</td><td>0.56</td><td>0.96</td><td>2.0</td></tr><tr><td>Radio Galaxies</td><td>3.4</td><td>4.8 × 10<sup>-3</sup></td><td>0.10</td><td>4.2</td><td>5.1</td><td>6.9 × 10<sup>-3</sup></td><td>0.13</td><td>4.7</td></tr><tr><td>Star-forming Galaxies</td><td>0.030</td><td>0.37</td><td>0.93</td><td>2.0</td><td>~0</td><td>~1</td><td>1.0</td><td>1.7</td></tr><tr><td>Dust-obscured AGNs</td><td>1.0 × 10<sup>-3</sup></td><td>0.73</td><td>0.98</td><td>1.5</td><td>~0</td><td>~1</td><td>1.0</td><td>1.4</td></tr><tr><td>IceCube High-energy Tracks</td><td>0.77</td><td>0.05</td><td>0.49</td><td>5.2</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** The test statistics values are reported as $\lambda$, the pre-trial $p$-values are labelled as $p$, and the post-trial values as $P$. The 90% CL upper limits on $\nu_{\mu} + \bar{\nu}_{\mu}$ flux are expressed in terms of the total $E^{-2}$ flux normalization at 1 GeV for the considered objects (in units of $10^{-8} \text{ GeV}^{-1} \text{ cm}^{-2} \text{ s}^{-1}$).

## 5. Results

The results of the stacking analysis are summarized in Table 2. The computation of the trial factor has been performed by generating $10^4$ pseudo-experiments with only background events (randomized samples), and applying the same likelihood analysis for the 10 different combinations of catalogs and weighting schemes considered in this analysis.

None of the tested catalogs shows a significant association with the 11 yr point-source sample of neutrinos detected by ANTARES. The smallest $p$-value is obtained for the radio galaxy catalog with the equal-weights hypothesis, with a pre-trial $p$-value $p = 4.8 \times 10^{-3}$, equivalent to a $2.8\sigma$ excess (two-sided convention), reduced to $1.6\sigma$ post-trial.

The 90% confidence level (CL) upper limits on flux in Table 2 are expressed in terms of the one-flavor ($\nu_{\mu} + \bar{\nu}_{\mu}$) total flux normalization at 1 GeV:

$$ \Phi_{\text{tot}}^{0} = \sum_{j=1}^{N_{\text{sources}}} \Phi_{j}^{0}, \qquad (10) $$

computed for a spectral index $\gamma = 2$. For comparison, in Huber (2019) the IceCube Collaboration reports the results of a similar stacking likelihood analysis using 1301 blazars contained in the Fermi 3FHL catalog (Ajello et al. 2017). The quoted upper limit for an $E^{-2}$ spectrum is $\Phi_{90\%}^{UL} = 7.2 \times 10^{-9} \text{ GeV}^{-1} \text{ cm}^{-2} \text{ s}^{-1}$, which is a factor ~6 more stringent than the ANTARES limit presented here.

### 5.1. *Search for Dominant Sources*

To search for a possible dominant contribution of a small number of sources, an individual likelihood fit has been performed for each object in the radio galaxy and Fermi 3LAC catalogs. The resulting $p$-values associated with each source are shown in Figure 4. Apart from the bulk of the distribution, a single source is present in the tail of the distribution of $p$-values for both the radio galaxies and the Fermi blazars. Their properties are thus investigated in more details in the following sections.

**Radio galaxies.** The radio galaxy with the smallest $p$-value is 3C 403, which is found with a pre-trial $p$-value of $p = 2.3 \times 10^{-4}$, equivalent to $3.7\sigma$. The probability of finding by chance a similar excess among the $N = 56$ radio galaxies in the field of view of ANTARES has been computed using $P(n \ge 1) = 1 - (1 - p)^N$, leading to a value $P = 1.3 \times 10^{-2}$, equivalent to $2.5\sigma$.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Distribution of the individual $p$-values for the Fermi 3LAC blazars and for the radio galaxies.</figcaption></figure>

The arrival directions of ANTARES events around this source are shown in Figure 5(left), where one can see that the excess comes from the presence of two events lying at less than $0.5^\circ$ from the source. Those track events have an estimated angular uncertainty of $0.2^\circ$, and reconstructed energies of ~5 TeV and ~10 TeV, respectively.

The radio galaxy is located at a declination of $+2.5^\circ$ and its distance is evaluated to be ~260 Mpc. The galaxy 3C 403 is of type FR-II and presents a peculiar radio morphology, called “X-shape” or “Winged” (see Figure 1 of Dennett-Thorpe et al. 2002) due to a recent change in the orientation of the jet axis, which could be the signature of a binary black hole merger or fusion with a small galaxy. The interaction of the newly formed jet with the radiation and matter environment around the black hole could lead to high-energy neutrino emission, as suggested in Kun et al. (2018).

**Blazars.** The most significant blazar from the Fermi 3LAC catalog is the BL Lac MG3 J225517+2409 (Fermi source 3FGL J2255.1+2411, or 4FGL J2255.2+2411 in the latest version of the catalog). A recent compilation of multi-wavelength observations for this object is reported in Figure 8 of Giommi et al. (2020), showing a synchrotron peak around $10^{15} \text{ Hz}$, indicating that the blazar can be classified in the ISP category, similar to TXS 0506+056. The redshift of MG3 J225517+2409 is currently unknown, but a limit of $z > 0.86$ has been recently reported in Paiano et al. (2019).

The gamma-ray flux of MG3 J225517+2409 measured by Fermi-LAT in eight years of observation (Abdollahi et al. 2020) is well described by a simple power law without cutoff,
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Arrival directions of ANTARES track events in equatorial coordinates around the most prominent sources: radio galaxy 3C 403 (left) and blazar MG3 J225517+2409 (right). Each event is represented by a dashed empty circle with a radius equal to its angular error estimate. The estimated energy of each event is indicated by a color code. The thin black circles centered around the source (red star marker) indicate the 1° and 5° angular distances to the source. On the right plot, the IceCube track IC-100608A is indicated as a green dashed circle, with radius equal to the 50% containment value (~1.2°).</figcaption></figure>

with a spectral index $\gamma = 2.10 \pm 0.05$ and an integral flux of $\Phi = 1.0 \times 10^{-9}$ photons cm$^{-2}$ s$^{-1}$ for energies between 1 and 100 GeV.

It is noteworthy that another candidate blazar, 3FGL J2258.8+2437, is located only ~1° away from MG3 J225517+2409, and ~1.2° away from the position of the IceCube high-energy track event IC-100608A. However, this source does not have an identified counterpart and is not present in the 10 yr Fermi-LAT catalogs, suggesting that 3FGL J2258.8+2437 might have been a spurious detection. For the present study, this source is thus not considered in the calculation of chance coincidence.

The result of the individual likelihood fit for this source gives a pre-trial *p*-value of $p = 1.4 \times 10^{-4}$, equivalent to $3.8\sigma$. The probability of finding by chance a similar excess among the $N = 1255$ blazars in the field of view of ANTARES is evaluated by generating pseudo-experiments, giving a post-trial value $P = 0.15$ (equivalent to $1.4\sigma$). The simple binomial evaluation gives a very similar value, $P(n \ge 1) = 1 - (1-p)^N \approx 0.16$, because the overlap between sources is very small: there are only 58 pairs of 3LAC objects that are separated by an angular distance less than 1°.

The distribution of ANTARES events around the blazar is shown in Figure 5(right). There are five ANTARES tracks lying at less than 1° from the source, with estimated energies ranging from ~3 to ~40 TeV and with estimated angular uncertainty between 0.2° and 0.5°. This particular region of the sky has been previously reported by the ANTARES Collaboration (Illuminati et al. 2019) to contain the hottest spot of the full-sky point-source blind search; the maximum excess lies at <0.7° from the position of the blazar MG3 J225517+2409.

Under the (restrictive) assumption that a neutrino flux is produced steadily by the source with an $E^{-2}$ spectrum, an estimation of the flux normalization at 1 GeV of MG3 J225517+2409 with the ANTARES data would be $\Phi_0 \sim (0.6-3.0) \times 10^{-8}$ GeV$^{-1}$ cm$^{-2}$ s$^{-1}$, depending on the number of neutrino candidates that are considered as signal (between 1 and 5) from this source.

In addition to the ANTARES measurement, a high-energy muon track IC-100608A from the IceCube 8 yr up-going muon sample (Haack & Wiebusch 2017) is found to lie at an angular distance of 1.1° from the source. This event has a deposited energy in the IceCube detector of about ~300 TeV, but it has a poor angular resolution: the 50% CL containment region has a radius ~1°, but the estimated 90% CL region has a radius larger than ~3°. This poor angular reconstruction explains why the spatial correlation with MG3 J225517+2409 has not been reported before, the IceCube muon event IC-100608A being excluded from the selection in previous blazar correlation searches (Padovani et al. 2016; Giommi et al. 2020).

## 5.2. Dedicated Analysis of the Region around the Blazar MG3 J225517+2409

The possible association between ANTARES and IceCube neutrinos and the blazar MG3 J225517+2409 is investigated in more detail in this section, by further looking at the time evolution of the gamma-ray emission of the source. The latest data from the Fermi 4FGL catalog (Abdollahi et al. 2020) are used to compare the evolution of the gamma-ray flux with the arrival times of ANTARES and IceCube events.

The Fermi 4FGL catalog provides measurements with twice as much exposure as the Fermi 3FGL, with 20% larger acceptance, improved angular resolution above 3 GeV, a better model of the background galactic diffuse emission, and a refined likelihood analysis. The 4FGL catalog was not available at the time of the stacking analysis; it is used only in this section for the time analysis.

The time evolution of the gamma-ray flux measured by Fermi-LAT for the blazar 4FGL J2255.2+2411 above 100 MeV is shown in Figure 6, along with the arrival times of spatially coincident IceCube and ANTARES neutrino events. The data are presented in 48 bins, each two months wide, from 2008 August to 2016 August. A Bayesian block algorithm (Scargle et al. 2013) has been applied to the flux data in order to pinpoint possible flaring states. A significant increase of the flux by a factor of ~4 with respect to the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Gamma-ray flux for the BL Lac 4FGL J2255.2+2411 measured by Fermi-LAT as a function of time in modified Julian date (black points). The red line represents the model gamma-ray light curve that is used in the time-dependent likelihood (obtained via a Bayesian block algorithm). The blue vertical lines indicate the arrival times of ANTARES events located at an angular distance smaller than 2° from the blazar. The vertical green line corresponds to the arrival time of the IceCube track IC-100608A.</figcaption></figure>

average value is evident in a time window of ~4 months, centered at about MJD = 55,400.

### 5.2.1. Antares Events

The arrival times of the five ANTARES events around the blazar MG3 J225517+2409 are clustered at earlier observation times, in a large time window encompassing the flare but with none of them actually occurring during the period of high gamma-ray activity (the two latest events were observed ~2 and ~3 months after the end of the gamma-ray flare). Therefore, including a time-dependent term in the likelihood analysis does not strengthen the potential association between ANTARES candidate neutrino events and the gamma-ray emission of the blazar.

Considering only the time information, a very simple and model-independent test is to compare the average time to the flare center ($t_F$)

$$ \tau = \frac{1}{N} \sum_{i=1}^{N=5} |t_i - t_F|, \qquad (11) $$

observed for the five ANTARES events to the distribution that is expected from random pseudo-experiments (following the true time distribution of ANTARES events).

The value obtained for the data is $\tau \approx 400$ days, and the fraction of pseudo-experiments having a lower or equal value is found to be $\approx 2\%$, thus leading to an estimated $p$-value of $p = 0.02$ ($2.3\sigma$).

### 5.2.2. IceCube Event

To evaluate the degree of space and time correlation between the gamma-ray flux of MG3 J225517+2409 and the neutrino candidate of IceCube, an additional likelihood analysis has been performed for this specific event.

The likelihood of Equation (1) is modified as follows:

$$ \ln \mathcal{L}(x|H_1) = \ln [\mu_s S(x_i) f_S(t_i) + \mu_b B(x_i) f_B(t_i)] - \mu_s - \mu_b, \qquad (12) $$

where the new terms $f_S(t_i)$ and $f_B(t_i)$ represent the signal and background time PDFs respectively.

The signal $f_S(t)$ is chosen to be proportional to the Bayesian block estimation (Scargle et al. 2013) of the gamma-ray flux of MG3 J225517+2409, shown as a solid red curve in Figure 6. The background $f_B(t)$ is proportional to the average number of events detected per day by IceCube, using the point-source selection cuts (Aartsen et al. 2019). The value of $f_B(t)$ is assumed to be constant during the different data-taking periods, from the IC59 sample (2009–2010) where $f_B(t) \approx 60$ events/day, up to the IC86 sample where $f_B(t) \approx 210$ events/day.

The spatial PDF $S(x)$ is taken as a two-dimensional Gaussian function of width $\sigma = 1.8^\circ$, corresponding to the 90% CL statistical angular error of ~3° reported in Haack & Wiebusch (2017). The background distribution $B(x)$ is assumed to be uniform in right ascension, and proportional to the value of a Bayesian block fit of the distribution in declination of the 56 IceCube high-energy tracks considered in the stacking analysis, as shown in Figure 7.

The result of the likelihood fit gives a fitted number of signal events $n_s = 0.68$ and a test statistic equal to $\lambda = 0.39$. The associated $p$-value is computed by generating pseudo-experiments, leading to a value $p = 5.3 \times 10^{-2}$, equivalent to $1.9\sigma$. As the considered source has been specified a priori using the ANTARES result, there are no trial factors to account for.

## 6. Discussion and Conclusion

The results have been presented of a stacking point-source likelihood search using 11 years of ANTARES data with muon neutrino candidates reconstructed as tracks. After accounting for trial factors, none of the catalogs considered have shown a significant result in the stacking analysis.

The most significant value is obtained for the radio galaxy catalog, with a pre-trial value of $p = 4.8 \times 10^{-3}$, reducing to $P = 0.1$ post-trial. Of the individual sources, the radio galaxy 3C 403 ($P = 0.013$) and the blazar MG3 J225517+2409 ($P = 0.15$) are reported as interesting targets.
<figure><img src="image_9.png" /><figcaption>Figure 7. Distribution of the sine of the declination of the 56 high-energy track-like events detected by IceCube considered in this study. The blue filled histogram represents the Bayesian block fit to the data that is used in the likelihood analysis.</figcaption></figure>

For the Fermi 3LAC blazars, the flux upper limits obtained in this study for an $E^{-2}$ neutrino spectrum and equal weights are a factor ~6 less stringent than the comparable IceCube limits presented in Huber (2019). This also limits the maximum contribution of blazars to the diffuse astrophysical flux in IceCube. Following the same procedure as Huber (2019), which quotes a maximum contribution of 13.0%–16.7%, our less-stringent limits yield a maximum contribution of ~80%–100% at 90% CL. Since ANTARES probes a lower energy range than that in which the IceCube flux is detected, we expect our limit to improve for steeper spectral indices, which are currently favoured by model fits (Aartsen et al. 2020b).

For the radio galaxies, the mild excess found in ANTARES data gives flux upper limits of the same order of magnitude as for blazars, allowing for a non-negligible contribution from this class of source to the diffuse astrophysical flux in IceCube.

In addition, a high-energy IceCube track is observed to be in coincidence with the blazar MG3 J225517+2409. This motivated an a posteriori analysis of the space and time correlation between the gamma-ray flux measured by Fermi and the neutrinos detected by ANTARES and IceCube for this source. The IceCube event is observed to be temporally coincident with what could be called a flaring period of the blazar. The chance probability of this association is estimated to be $p = 5.3 \times 10^{-2}$ (1.9σ).

The five ANTARES events lying within 1° from the source are not coincident in time with the flare. However, they all arrive within 400 days of the flare, with an associated p-value of $p = 0.02$. This result hardly reconciles a direct association between a TeV neutrino excess and an MeV–GeV gamma-ray flare. A similar finding was already highlighted for the apparent lack of association between the 2014/15 neutrino excess observed by IceCube from TXS 0506+056 and observation of Fermi-LAT activity (Aartsen et al. 2018b; Garrappa & Buson 2019).

The lack of gamma rays in coincidence with the ANTARES neutrino events could be due to a high optical depth to photons at Fermi-LAT energies. Given the current range of models for neutrino production in TXS 0506+056 (Keivani et al. 2018; Reimer et al. 2019; Cao et al. 2020; Sahu et al. 2020; Petropoulou et al. 2020; Murase et al. 2018), such a scenario seems plausible.

Combining the space and time probabilities for the ANTARES events using

$$p_{\text{space-time}} = p_{\text{spatial}} \times p_{\text{time}} (1 - \ln(p_{\text{spatial}} \times p_{\text{time}})),$$

leads to a value $p_{\text{space-time}} = 1.5 \times 10^{-2}$ (2.3σ).

Combining the ANTARES spacetime and the IceCube p-values gives

$$\begin{aligned} p_{\text{combined}} &= p_{\text{ANT}} \times p_{\text{IC}} (1 - \ln(p_{\text{ANT}} \times p_{\text{IC}})) \\ &= 8.0 \times 10^{-3} (2.6\sigma). \end{aligned}$$

Thus there is no strong statistical evidence for an association. Furthermore, we emphasize that these tests are done a posteriori, and that the quoted p-values reported here are thus subject to bias. Nonetheless, given the suggestions of neutrino associations with other blazars (Kadler et al. 2016; Garrappa & Buson 2019), we advocate to include MG3 J225517+2409 as a potential interesting neutrino source for future studies.

The accumulation of independent neutrino data with larger statistics and better energy and angle resolution in the near future by next-generation neutrino telescopes such as KM3NeT (Adrián-Martínez et al. 2016) will be of great interest for studying the correlation between gamma-ray and neutrino emission.

The authors thank Loredana Bassani for the suggestion of the radio galaxy catalog and for fruitful discussions.

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), Institut Universitaire de France (IUF), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Région Île-de-France (DIM-ACAV), 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; Council of the President of the Russian Federation for young scientists and leading scientific schools supporting grants, Russia; Executive Unit for Financing Higher Education, Research, Development and Innovation (UEFISCDI), Romania; Ministerio de Ciencia, Innovación, Investigación y Universidades (MCIU): Programa Estatal de Generación de Conocimiento (refs. PGC2018-096663-B-C41, -A-C42, -B-C43, -B-C44) (MCIU/FEDER), Severo Ochoa Centre of Excellence and MultiDark Consolider (MCIU), Junta de Andalucía (ref. SOMM17/6104/UGR and A-FQM-053-UGR18), Generalitat Valenciana: Grisolía (ref. GRISOLIA/2018/119) and GenT (ref. CIDEGENT/2018/034), Spain; Ministry of Higher Education, Scientific Research and Professional Training, Morocco. We also acknowledge the technical support of Ifremer, AIM and Foselev Marine for the sea operation and the CC-IN2P3 for the computing facilities.

## ORCID iDs

J. Aublin https://orcid.org/0000-0001-7642-6225
B. Baret https://orcid.org/0000-0001-6064-3858
S. Basa https://orcid.org/0000-0002-4291-333X
S. Biagi https://orcid.org/0000-0001-8598-0017
M. Bissinger https://orcid.org/0000-0002-8709-8236

S. Celli https://orcid.org/0000-0002-7592-0851

L. Fusco https://orcid.org/0000-0001-8254-3372

J. Hofestädt https://orcid.org/0000-0002-7848-117X

G. Illuminati https://orcid.org/0000-0002-4138-8027

C. W. James https://orcid.org/0000-0002-6437-6176

M. Kadler https://orcid.org/0000-0001-5606-6154

O. Kalekin https://orcid.org/0000-0001-6206-1288

A. Kouchner https://orcid.org/0000-0001-7068-2113

I. Kreykenbohm https://orcid.org/0000-0001-7335-1803

E. Leonora https://orcid.org/0000-0002-0536-3551

A. Margiotta https://orcid.org/0000-0001-6929-5386

A. Marinelli https://orcid.org/0000-0002-1466-1219

S. Navas https://orcid.org/0000-0003-1688-5758

P. Piattelli https://orcid.org/0000-0003-4748-6485

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

S. J. Tingay https://orcid.org/0000-0002-8195-7562

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

J. Wilms https://orcid.org/0000-0003-2065-5410

## References

Aab, A., Abreu, P., Aglietta, M., et al. 2018, ApJ, 853, L29

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, ApJ, 835, 45

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

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

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

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, EPJC, 79, 234

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020a, PhRvL, 124, 051103

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

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Ackermann, M., Ajello, M., Allafort, A., et al. 2011, ApJ, 743, 171

Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 755, 164

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015, ApJ, 810, 14

Adrián-Martínez, S., Ageron, M., Aharonian, F., et al. 2016, JPhG, 43, 084001

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

Ajello, M., Atwood, W. B., Baldini, L., et al. 2017, ApJS, 232, 18

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

Albert, A., André, M., Anghinolfi, M., et al. 2018, ApJL, 863, L30

Albert, A., André, M., Anghinolfi, M., et al. 2020, ApJ, 892, 92

Barlow, R. 1990, NIMPA, 297, 496

Bassani, L., Venturi, T., Molina, M., et al. 2016, MNRAS, 461, 3165

Cao, G., Yang, C., Yang, J., et al. 2020, PASJ, 72, 20

Dennett-Thorpe, J., Scheuer, P., Laing, R. A., et al. 2002, MNRAS, 330, 609

Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433

Garrappa, S., Buson, S., Franckowiak, A., et al. 2019, ApJ., 880, 103

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

Haack, C., & Wiebusch, C. 2017, ICRC (Busan), 35, 1005

Halzen, F., & Zas, E. 1997, ApJ, 488, 669

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hooper, D. 2016, JCAP, 2016, 002

Huber, M. 2019, ICRC (Madison, WI), 36, 916

Huber, M., & Krings, K. 2017, ICRC (Busan), 35, 994

Illuminati, G., Aublin, J., Navas, S., et al. 2019, ICRC (Madison, WI), 36, 920

Kadler, M., Krauß, F., Mannheim, K., et al. 2016, NatPh, 12, 807

Keivani, A., Murase, K., Petropoulou, M., et al. 2018, ApJ, 864, 84

Kopper, C. 2017, ICRC (Busan), 35, 981

Kun, E., Biermann, P., Britzen, S., et al. 2018, Univ, 4, 24

Maggi, G., Buitink, S., Correa, P., et al. 2016, PhRvD, 94, 103007

Mannheim, K. 1995, APh, 3, 295

Mücke, A., Protheroe, R. J., Engel, R., et al. 2003, APh, 18, 593

Murase, K., Inoue, Y., & Dermer, C. D. 2014, PhRD, 90, 023007

Murase, K., Oikonomou, F., & Petropoulou, M. 2018, ApJ, 865, 124

Padovani, P., & Giommi, P. 1995, ApJ, 444, 567

Padovani, P., Resconi, E., Giommi, P., et al. 2016, MNRAS, 457, 3582

Paiano, S., Padovani, P., Falomo, R., et al. 2019, ATel, 13202, 1

Petropoulou, M., Murase, K., Santande, M., et al. 2020, ApJ, 891, 115

Reimer, A., Böttcher, M., & Buson, S. 2019, ApJ, 881, 46

Sahu, S., López Fortín, C. E., & Nagataki, S. 2020, ApJ, 898, 103

Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167

Schnabel, J. 2013, NIMPA, 725, 106

Tavecchio, F., & Ghisellini, G. 2015, MNRAS, 451, 1502

van Velzen, S., Falcke, H., Schellart, P., et al. 2012, A&A, 544, A18