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

CrossMark

# A Search for Cosmic Neutrino and Gamma-Ray Emitting Transients in 7.3yr of ANTARES and Fermi LAT Data

H. A. Ayala Solares¹,², D. F. Cowen¹,²,³, J. J. DeLaunay¹,² ID, D. B. Fox²,³,⁴ ID, A. Keivani¹,² ID, M. Mostafá¹,²,³, K. Murase¹,²,³ ID, C. F. Turley¹,² ID

## AMON,

and

A. Albert⁵ ID, M. André⁶, M. Anghinolfi⁷, G. Anton⁸, M. Ardid⁹, J.-J. Aubert¹⁰, J. Aublin¹¹, B. Baret¹¹ ID, J. Barrios-Martí¹², S. Basa¹³ ID, B. Belhorma¹⁴, V. Bertin¹⁰, S. Biagi¹⁵, R. Bormuth¹⁶,¹⁷, J. Boumaaza¹⁸, S. Bourret¹⁹, M. Bouta²⁰, M. C. Bouwhuis¹⁶, H. Brânzaş²¹, R. Bruijn¹⁶,²², J. Brunner¹⁰, J. Busto¹⁰, A. Capone²³,²⁴, L. Caramete²¹, J. Carr¹⁰, S. Celli²³,²⁴,²⁵, M. Chabab²⁶, R. Cherkaoui El Moursli¹⁸, T. Chiarusi²⁷, M. Circella²⁸, A. Coleiro¹¹,¹², M. Colomer¹¹,¹², R. Coniglione¹⁵, H. Costantini¹⁰, P. Coyle¹⁰, A. Creusot¹¹, A. F. Díaz²⁹, A. Deschamps³⁰, C. Distefano¹⁵, I. Di Palma²³,²⁴, A. Domi⁷,³¹, R. Donà²⁷, C. Donzaud¹¹,³², D. Dornic¹⁰, D. Drouhin⁵, T. Eber¹⁸, I. El Bojaddaini²⁰, N. El Khayati¹⁸, D. Elsässer³³, A. Enzenhöfer⁸,¹⁰, A. Ettahiri¹⁸, F. Fassi¹⁸, P. Fermani²³,²⁴, G. Ferrara¹⁵, L. Fusco¹¹,³⁴ ID, P. Gay¹¹,³⁵, H. Glotin³⁶, R. Gozzini¹², T. Grégoire¹¹, R. Gracia Ruiz⁵, K. Graf⁸, S. Hallmann⁸, H. van Haren³⁷, A. J. Heijboer¹⁶, Y. Hello³⁰, J. J. Hernández-Rey¹², J. Höfl⁸, J. Hofestädt⁸, G. Illuminati¹² ID, C. W. James³⁸,³⁹, M. de Jong¹⁶,¹⁷, M. Jongen¹⁶, M. Kadler³³ ID, O. Kalekin⁸ ID, U. Katz⁸ ID, N. R. Khan-Chowdhury¹², A. Kouchner¹¹,⁴⁰, M. Kreter³³, I. Kreykenbohm⁴¹, V. Kulikovskiy⁷,⁴², R. Lahmann⁸, R. Le Breton¹¹ ID, D. Lefèvre⁴³,⁴⁴, E. Leonora⁴⁵ ID, G. Levi²⁷,³⁴, M. Lincetto¹⁰, D. Lopez-Coto⁴⁶, M. Lotze¹², S. Loucatos¹¹,⁴⁷, G. Maggi¹⁰, M. Marcelin¹³, A. Margiotta²⁷,³⁴, A. Marinelli⁴⁸,⁴⁹ ID, J. A. Martínez-Mora⁹, R. Mele⁵⁰,⁵¹ ID, K. Melis¹⁶,²², P. Migliozzi⁵⁰, A. Moussa²⁰, S. Navas⁴⁶, E. Nezri¹³, C. Nielsen¹¹, A. Nuñez¹⁰,¹³, M. Organokov⁵, G. E. Päväläns²¹, C. Pellegrino²⁷,³⁴, M. Perrin-Terrin¹⁰, P. Piattelli¹⁵, V. Popa²¹, T. Pradier⁵, L. Quinn¹⁰, C. Racca⁵², N. Randazzo⁴⁵, G. Riccobene¹⁵, A. Sánchez-Losa²⁸, A. Salah-Eddine²⁶, I. Salvadori¹⁰, D. F. E. Samtleben¹⁶,¹⁷, M. Sanguineti⁷,³¹ ID, P. Sapienza¹⁵, F. Schüssler⁴⁷ ID, M. Spurio²⁷,³⁴ ID, Th. Stolarczyk⁴⁷ ID, M. Taiuti⁷,³¹, Y. Tayalati¹⁸, T. Thakore¹², A. Trovato¹⁵, B. Vallage¹¹,⁴⁷, V. Van Elewyck¹¹,⁴⁰, F. Versari²⁷,³⁴, S. Viola¹⁵, D. Vivolo⁵⁰,⁵¹, J. Wilms⁴¹ ID, D. Zaborov¹⁰, J. D. Zornoza¹², and J. Zúñiga¹²

## ANTARES Collaboration

¹Department of Physics, Pennsylvania State University, University Park, PA 16802, USA; cft114@psu.edu

²Center for Particle & Gravitational Astrophysics, Institute for Gravitation and the Cosmos, Pennsylvania State University, University Park, PA 16802, USA

³Department of Astronomy & Astrophysics, Pennsylvania State University, University Park, PA 16802, USA

⁴Center for Theoretical & Observational Cosmology, Institute for Gravitation and the Cosmos, Pennsylvania State University, University Park, PA 16802, USA

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

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

⁷INFN-Sezione di Genova, Via Dodecaneso 33, I-16146 Genova, Italy

⁸Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, D-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, E-46730 Gandia, Spain

¹¹APC, Univ Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs de Paris, Sorbonne Paris Cité, France

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

¹³LAM—Laboratoire d'Astrophysique de Marseille, Pôle de l'Étoile Site de Château-Gombert, rue Frédéric Joliot-Curie 38, F-13388 Marseille Cedex 13, 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, I-95123 Catania, Italy

¹⁶Nikhef, Science Park, Amsterdam, The Netherlands

¹⁷Huygens-Kamerlingh Onnes Laboratorium, Universiteit Leiden, The Netherlands

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

¹⁹PC, Univ Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs de Paris, Sorbonne Paris Cité, France

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

²¹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, I-00185 Roma, Italy

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

²⁵Gran Sasso Science Institute, Viale Francesco Crispi 7, I-00167 L'Aquila, Italy

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

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

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

²⁹Department of Computer Architecture and Technology/CITIC, University of Granada, E-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, I-16146 Genova, Italy

³²Université Paris-Sud, F-91405 Orsay Cedex, France

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

³⁴Dipartimento di Fisica e Astronomia dell'Università, Viale Berti Pichat 6/2, I-40127 Bologna, 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, F-83041 Toulon, France

³⁷Royal Netherlands Institute for Sea Research (NIOZ) and Utrecht University, Landsdip 4, 1797 SZ 't Horntje (Texel), The Netherlands

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

³⁹ARC Centre of Excellence for All-sky Astrophysics (CAASTRO), Australia
⁴⁰ Institut Universitaire de France, F-75005 Paris, France

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

⁴² Moscow State University, Skobeltsyn Institute of Nuclear Physics, Leninskie gory, 119991 Moscow, Russia

⁴³ Mediterranean Institute of Oceanography (MIO), Aix-Marseille University, F-13288, Marseille, Cedex 9, France

⁴⁴ Université du Sud Toulon-Var, CNRS-INSU/IRD UM 110, F-83957, La Garde Cedex, France

⁴⁵ INFN—Sezione di Catania, Via S. Sofia 64, I-95123 Catania, Italy

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

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

⁴⁸ INFN—Sezione di Pisa, Largo B. Pontecorvo 3, I-56127 Pisa, Italy

⁴⁹ Dipartimento di Fisica dell'Università, Largo B. Pontecorvo 3, I-56127 Pisa, Italy

⁵⁰ INFN—Sezione di Napoli, Via Cintia I-80126 Napoli, Italy

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

⁵² GRPHE—Université de Haute Alsace—Institut universitaire de technologie de Colmar, 34 rue du Grillenbreit BP, F-50568-68008 Colmar, France
Received 2019 April 12; revised 2019 September 30; accepted 2019 October 1; published 2019 November 25

## Abstract

We analyze 7.3 yr of ANTARES high-energy neutrino and Fermi Large Area Telescope (LAT) γ-ray data in search of cosmic neutrino + γ-ray (ν+γ) transient sources or source populations. Our analysis has the potential to detect either individual ν+γ transient sources (durations δt ≲ 1000 s), if they exhibit sufficient γ-ray or neutrino multiplicity, or a statistical excess of ν+γ transients of individually lower multiplicities. Individual high γ-ray multiplicity events could be produced, for example, by a single ANTARES neutrino in coincidence with a LAT-detected γ-ray burst. Treating ANTARES track and cascade event types separately, we establish detection thresholds by Monte Carlo scrambling of the neutrino data, and determine our analysis sensitivity by signal injection against these scrambled data sets. We find our analysis is sensitive to ν+γ transient populations responsible for >5% of the observed gamma-coincident neutrinos in the track data at 90% confidence. Applying our analysis to the unscrambled data reveals no individual ν+γ events of high significance; two ANTARES track + Fermi γ-ray events are identified that exceed a once per decade false alarm rate threshold (p = 17%). No evidence for subthreshold ν+γ source populations is found among the track (p = 39%) or cascade (p = 60%) events. Exploring a possible correlation of high-energy neutrino directions with Fermi γ-ray sky brightness identified in previous work yields no added support for this correlation. While TXS 0506+056, a blazar and variable (nontransient) Fermi γ-ray source, has recently been identified as the first source of high-energy neutrinos, the challenges in reconciling observations of the Fermi γ-ray sky, the IceCube high-energy cosmic neutrinos, and ultrahigh-energy cosmic rays using only blazars suggest a significant contribution by other source populations. Searches for transient sources of high-energy neutrinos thus remain interesting, with the potential for either neutrino clustering or multimessenger coincidence searches to lead to discovery of the first ν+γ transients.

**Key words:** BL Lacertae objects: general – cosmic rays – gamma-ray burst: general – gamma rays: general – neutrinos

## 1. Introduction

The ANTARES telescope (Ageron et al. 2011) is a deep-sea Cerenkov neutrino detector, located 40 km offshore from Toulon, France, in the Mediterranean Sea. The detector comprises a three-dimensional array of 885 optical modules, each one housing a 10 inch photomultiplier tube, and distributed over 12 vertical strings anchored in the sea floor at a depth of about 2400 m. The detection of light from upgoing charged particles is optimized with the photomultipliers facing 45° downward. Completed in 2008 May, the telescope aims primarily at the detection of neutrino-induced muons that cause the emission of Cerenkov light in the detector (track-like events). Charged current interactions induced by electron neutrinos (and, possibly, by tau neutrinos of cosmic origin) or neutral current interactions of all neutrino flavors can be reconstructed as cascade-like events (Albert et al. 2017a).

Due to its location, the ANTARES detector mainly observes the Southern sky (2π sr at any time). Events arising from sky positions in the decl. band −90° ≤ δ ≤ −48° are always visible as upgoing. Neutrino-induced events in the decl. band −48° ≤ δ ≤ +48° are visible as upgoing with a fraction of time decreasing from 100% down to 0%. While ANTARES has a substantially smaller volume than IceCube, the use of sea water as detection medium (rather than ice) provides better pointing resolution for individual events, especially those of cascade type, and its geographic location enables reduced-

background studies of the Southern hemisphere including the Galactic center region. On the other hand, natural light emission in the water leads to higher background levels (ANTARES Collaboration et al. 2005).

Chief scientific results from ANTARES include: searches for neutrino sources using track- and cascade-like events in data collected between 2007 and 2015 (Albert et al. 2017b); dedicated studies along the Galactic plane (Albert et al. 2017c), also in collaboration with the IceCube telescope (Albert et al. 2018a); and searches for an excess of high-energy cosmic neutrinos over the background of atmospheric events (Albert et al. 2018b). No cosmic neutrinos have been positively identified in the ANTARES data. Despite this, by integrating the cosmic neutrino spectrum from IceCube Collaboration et al. (2017) over the ANTARES effective area (Albert et al. 2017b), we estimate an expected 6.8 neutrinos of cosmic origin are detected each year, though all but the most energetic will be indistinguishable from the atmospheric background. Among all the possible astrophysical sources, transient sources increase the observation possibilities thanks to the suppression of atmospheric background in a well-defined spacetime window. For this reason, the ANTARES Collaboration is involved in a broad multimessenger program to exploit the connection between neutrinos and other cosmic messengers, including follow-up analyses associated with gravitational wave events
(Albert et al. 2017d; Albert et al. 2019b); coincidence searches against electromagnetic observations from radio (Croft et al. 2016; Albert et al. 2019a) and visible (Adrián-Martínez et al. 2016) to X- and γ-rays (Ageron et al. 2012); blazar flare episodes (Adrian-Martínez et al. 2015); and the neutrino source TXS 0506+056 (Albert et al. 2018c). To date, there have been no high-confidence counterparts identified for any ANTARES neutrino event.

In parallel, members of the Astrophysical Multimessenger Observatory Network (AMON;<fnref n="53" /> Smith et al. 2013; Cowen et al. 2016) have been exploring the possibility of neutrino + γ-ray (ν+γ) source identification via coincidence analysis, publishing analyses of Fermi Large Area Telescope (LAT; Atwood et al. 2009) and public IceCube 40-string (Keivani et al. 2015) and 59-string (Turley et al. 2018) data. Although no high-confidence ν+γ transients, nor evidence of subthreshold ν+γ source populations, were identified in these works, the latter revealed mild evidence for correlation between IceCube neutrino positions and the Fermi γ-ray sky.

Within the last year, a coincidence between the neutrino IceCube-170922A (Kopper et al. 2017) and the flaring blazar TXS 0506+056 (Tanaka et al. 2017) led to multimessenger (IceCube Collaboration et al. 2018b) and time-dependent neutrino clustering (IceCube Collaboration et al. 2018a) analyses suggesting this BL Lac-type object as the first known source of high-energy neutrinos and the first identified extragalactic cosmic ray accelerator. Further blazar source identifications can certainly be anticipated; however, the absence of point-source excesses in the ANTARES (Albert et al. 2017b) and IceCube (Aartsen et al. 2017a; Albert et al. 2018a) time-integrated data sets set strict limits on the fraction of cosmic high-energy neutrinos that can originate in these observed sources.

Possible alternative source populations include star-forming galaxies; starburst galaxies; galaxy groups and clusters; supernovae; and standard and low-luminosity gamma-ray bursts (see Murase 2015 for a review). Of these source possibilities, the transient and highly variable source populations will likely require time-sensitive searches for identification. Hadronic models foresee that neutrinos and γ-rays are cogenerated through the production and subsequent decay of mesons, mainly pions. γ-rays then result from the decay of neutral pions, while the decay of charged pions produces neutrinos. Additional processes in dense astrophysical regions can then degrade the energy of individual γ-rays to lower energies while leaving the neutrino energy spectrum almost unaffected, resulting in correlated emission of higher-energy neutrinos and lower-energy γ-rays.

The present paper is organized as follows. Details of the data sets are provided in Section 2. Our statistical approach and signal injection studies are discussed in Section 3. Unscrambled results and interpretation are presented in Section 4, and our conclusions in Section 5.

## 2. Data Sets

The Fermi LAT data set is highly complementary for cross-reference with high-energy neutrino data sets. The LAT offers a 1.4 steradian field of view, provides all-sky coverage every three hours on average, and exhibits good sensitivity over the 100 MeV $\lesssim \varepsilon_\gamma \lesssim 300$ GeV energy band.

This analysis was performed using publicly available Fermi LAT data. The relevant Fermi data were the Pass 8 photon reconstructions available from the LAT FTP server.<fnref n="54" /> These photon events were filtered using the Fermi Science Tools, keeping only photons with a zenith angle smaller than 90°, energies between 100 MeV and 300 GeV, detected during good time intervals as provided in the LAT satellite files.<fnref n="55" />

The point-spread function (PSF) of the LAT is given by a so-called double King function (King 1962) with the parameters depending on the photon energy, conversion type, and incident angle with respect to the LAT boresight (Ackermann et al. 2013). At energies in the hundreds of MeV, the angular uncertainty can be several degrees, especially for off-axis photons. At $\varepsilon_\gamma > 1$ GeV the average uncertainty drops below 1°, and at $\varepsilon_\gamma \gtrsim 100$ GeV angular uncertainties are better than 0°.1.

The ANTARES data used spans from 2007 February to 2015 December. Data from this 8.9 yr interval are divided into track and cascade events, all of which are upgoing. According to the selection criteria defined in Albert et al. (2017b), during this period 7622 track and 180 cascade neutrino candidates were identified. The Fermi mission has public data available starting from 2008 August 4. The ANTARES data is coincident with weeks 9 through 396 of the Fermi data, with 6774 track-like events and 162 cascade-like events falling within that 7.3 yr window. For the ANTARES data, the average PSFs for tracks and cascades are derived from Monte Carlo simulation, and then interpolated. For track and cascade events, the 90% containment radii for the PSFs are 1°.5 and 10°, respectively.

A healpix (Górski et al. 2005) map of resolution 8 (NSide = 256, mean spacing of 0°.23) was constructed using the entire Fermi data set (weeks 9 to 495 at the time of creation) with aforementioned photon selection criteria. Using the HEASoft software,<fnref n="56" /> events were binned into three logarithmically uniform energy bins. Each energy bin was then further binned into a healpix map, with the live time calculated via a Monte Carlo simulation. Dividing the counts map by the live time map produced the Fermi exposure map. Zero-valued (low-exposure) pixels were replaced by the average of the nearest neighbor pixels. Our three resulting all-sky Fermi maps are shown in Figure 1. Due to the additional reconstruction uncertainty in the Fermi PSF for high-inclination events (inclination angle greater than 60°), three additional maps for analysis of these events were generated by further averaging all pixels with their nearest neighbors.

## 3. Methods

### 3.1. Significance Calculation

Our analysis follows as an extension to the methods presented in Turley et al. (2018). Different from previous work, our analysis allows for coincidences with both multiple photons and multiple neutrinos. Our analysis also covers both the track and cascade events detected by ANTARES. For track-like events, we use an angular acceptance window of 5°, while for cascade-like events, we use a 10° window. For both event types, the temporal acceptance window is $\pm 1000$ s. Neutrino multiplets are constrained to have each neutrino within both the

<footnote n="53">⁵³ AMON website: http://www.amon.psu.edu/.</footnote>

<footnote n="54">⁵⁴ LAT data located at ftp://legacy.gsfc.nasa.gov/fermi/data/lat/weekly/photon/.</footnote>

<footnote n="55">⁵⁵ Fermi satellite files located at ftp://legacy.gsfc.nasa.gov/fermi/data/lat/weekly/spacecraft/.</footnote>

<footnote n="56">⁵⁶ HEASoft website: https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Background maps of the Fermi LAT γ-ray sky. Fermi data are split into three logarithmically uniform bins in energy and divided by the mission-averaged exposure map for that energy range. Grayscale intensity encodes the resulting mission-averaged photon flux over each band in units of photons per 200 s m⁻² deg⁻².</figcaption></figure>

angular and temporal separation of each other neutrino. Photons must fall within the angular and temporal window as measured from the average neutrino position and time. For each coincidence, a pseudo-log-likelihood test statistic, λ, is calculated as follows:

$$ \lambda = 2 \ln \frac{P_{i\gamma}(\mathbf{x}) n_{\nu}! n_{\gamma}! \prod_{\nu,\gamma} \tau(\Delta t_i)}{\prod_{\gamma} B_{\gamma,i}(\mathbf{x})} + \sum_{\nu} \ln \frac{1 - p_{c,i}}{p_{c,i}}, \quad (1) $$

where $P_{i\gamma}$ is the product of the PSF of each LAT photon and each ANTARES neutrino at the best position, $\mathbf{x}$, with each PSF normalized to have units of probability per square degree. The LAT PSF for each photon additionally depends on the photon energy, inclination angle, and conversion type. In general, the closer the PSF centers are, the larger the resulting λ value. The $n_ν$ and $n_γ$ terms are, respectively, the number of neutrinos and γ-rays in the coincidence. The $\Pi_{\nu,\gamma} \tau(\Delta t_i)$ term is the product of the temporal weighting function (Figure 2) evaluated for each neutrino and γ-ray in the coincidence.

For particles within 100 s of the average arrival time, this function is identically one, while it scales as $1/\Delta t$ for times between 100 and 1000 s. This allows the search to address the possibility of longer timescale associations (as might result from low-luminosity gamma-ray bursts; GRBs) while maintaining a preference for shorter timescale associations, if and when they are also present.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Temporal weighting function $\tau(\Delta t)$ used in the analyses. For $|\Delta t| < 100$ s, the function is flat and equal to 1. For $100 s < |\Delta t| < 1000$ s, the function scales as $1/\Delta t$.</figcaption></figure>

The $\Pi_{\gamma} B_{\gamma,i}(\mathbf{x})$ term is the product of LAT γ-ray backgrounds for each photon at the coincidence location, taken from the background maps shown in Figure 1. Together with the factorial terms, this acts like a Poisson probability of observing $n_γ$ photons from background. The $p_c$ factor, similar to the IceCube signalness (Aartsen et al. 2017b), is an energy proxy calculated by the ANTARES collaboration. The $p_c$ for a neutrino event is computed on an event-by-event basis using the normalised anti-cumulative distribution of the number of hits from the full ANTARES 2012–2017 neutrino data set. This probability represents the fraction of ANTARES events with a number of hits larger than that observed for the event: the larger the number of hits, the smaller the $p_c$ value. Overall, larger values of the λ statistic suggest a greater likelihood of a physically associated multiplet from a cosmic source, rather than a coincidence of uncorrelated events.

The best-fit position $\mathbf{x}$ is numerically calculated as the location of maximum PSF overlap. The photon multiplicity of each coincidence is calculated iteratively: beginning with a coincidence including all photons passing the temporal and proximity cuts, the photon with the lowest PSF density at the best-fit position is removed and a new λ, for the new best-fit position, is calculated. This process is repeated until one photon is left ($n_γ$ iterations), with the iteration yielding the maximum λ selected as the coincidence multiplicity.

This analysis presents two ways to identify a potential signal. First, with λ unbounded, the null distribution provides threshold values which can be used to identify individually significant coincidences and calculate their estimated false alarm rates. In this work, we use two such thresholds, $\lambda_D$ and $\lambda_C$, corresponding to false alarm rates of one per decade and one per century, respectively. Second, the presence of a subthreshold population of $\nu+\gamma$ emitting sources can be identified by a difference in the cumulative distributions of λ values between the observed and scrambled (null) populations. By design, true coincidences will be biased to higher λ values, and a population containing a sufficient number of signal events can be distinguished from the null distribution via an Anderson-Darling k-sample test (Scholz & Stephens 1987).
<table><caption>**Table 1**

Coincidence Search Results</caption><thead><tr><th rowspan="2">Data Set</th><th rowspan="2">⟨n<sub>ν+γ</sub>⟩</th><th colspan="4">Thresholds</th><th colspan="3">Observed Values</th></tr><tr><th>λ<sub>D</sub></th><th>λ<sub>C</sub></th><th>n<sub>inj,1%</sub></th><th>n<sub>inj,0.1%</sub></th><th>n<sub>ν+γ</sub></th><th>λ<sub>max</sub></th><th>p<sub>A-D</sub></th></tr></thead><tbody><tr><td>Tracks, 100 s</td><td>2716 ± 36</td><td>18.5</td><td>25.4</td><td>205</td><td>260</td><td>2734</td><td>18.94</td><td>39%</td></tr><tr><td>1000 s</td><td>”</td><td>”</td><td>”</td><td>220</td><td>285</td><td>”</td><td>”</td><td>”</td></tr><tr><td>Cascades</td><td>83.6 ± 5.8</td><td>8.1</td><td>14.6</td><td>-</td><td>-</td><td>80</td><td>2.7</td><td>60%</td></tr><tr><td>Track Multiplets</td><td>0.48 ± 0.69</td><td>-</td><td>-9.3</td><td>-</td><td>-</td><td>0</td><td>-</td><td>-</td></tr></tbody></table>

**Note:** $⟨n_{ν+γ}⟩$ is the expected number of neutrinos observed in coincidence with one or more γ-rays, as derived from 10,000 Monte Carlo scrambled realizations of each data set. $λ_D$ and $λ_C$ are the thresholds above which a coincidence is observed only once per simulated decade or century, respectively. $n_{inj,1%}$ and $n_{inj,0.1%}$ are the number of injected signal events required in simulations to give Anderson–Darling test (Scholz & Stephens 1987) p-values of $p < 1\%$ and $p < 0.1\%$, respectively, by comparison to the null distributions for each data set. $n_{ν+γ}$ is the number of neutrinos observed in coincidence with one or more γ-rays in unscrambled data, $λ_{max}$ is the maximum observed $λ$ for each data set, and $p_{A-D}$ is the Anderson–Darling test p-value from comparison of the observed $λ$ distribution to the associated null distribution. Cells with a “-” could not be calculated, for reasons detailed in the main text.

## 3.2. Background Generation

We generate a set of 10,000 Monte Carlo scrambled versions of each of our data sets to characterize their null distributions and define analysis thresholds, prior to performing any study of the unscrambled data sets. Our scrambling procedure begins by first converting the coordinates of each neutrino to detector coordinates. The arrival time and azimuthal angle of each original neutrino $ν_i$ are then exchanged with another randomly selected neutrino $ν_j$. Each neutrino retains its original elevation. Finally, the coordinates are converted back to the equatorial system. This approach is similar to the method used in our previous work (Turley et al. 2016), with the primary difference being the use of detector coordinates for the scrambling procedure. Fermi LAT photons are not scrambled as the LAT data contains known sources and extensive (complex) structure. Coincidence analysis is carried out for each scrambled data set and $λ$ values are calculated for the resulting $ν+γ$ coincidences via Equation (1). Thresholds from this analysis for false alarm rates of 1 per decade ($λ_D$) and 1 per century ($λ_C$) are presented in Table 1.

In contrast to previous work (Turley et al. 2018), due to the sensitivity to multineutrino events and the use of both track and cascade events, we split the analysis into three separate parts. The first part is to detect all coincidences with single-neutrino track-like events. The second looks for coincidences with multineutrino track-like events. The third and final part is a search for coincidences with all single-neutrino cascade-like events. Multineutrino cascades are not considered, as there are no cascade-like events within the temporal acceptance window of each other.

## 3.3. Signal Injection

To estimate the sensitivity of our analysis to subthreshold populations of cosmic $ν+γ$ emitting sources, we generate a population of signal-like events. These events are injected into the scrambled data sets so that the injected distributions can be compared with the null distribution.

We determine the multiplicity of a generated signal event following the methods used in Turley et al. (2018). This method assumes a population of sources emitting one neutrino, with associated photon fluence distributed according to $N(S \ge S_0) \propto S_0^{-3/2}$. In this formulation, $N(S \ge S_0)$ is the number of events observed with a fluence greater than the threshold fluence $S_0$. Setting this minimum to 0.001 photons, we can invert this relationship and generate the expectation

value for the multiplicity of an arbitrary event in terms of a uniform random variable $u$ as $⟨n_γ⟩ = S_0 u^{-2/3}$. The distribution of $n_γ$ is then calculated by drawing randomly from a Poisson distribution with the expectation value $⟨n_γ⟩$. Excluding events with zero photons, this yields the following $n_γ$ distribution: 93.8% singlet, 4.5% doublet, 0.9% triplet, and 0.38%, 0.19%, 0.095%, 0.0567%, 0.0365%, 0.0244%, and 0.0174% for multiplicities four through 10.

A signal event of photon multiplicity $n_γ$ is then generated by choosing a random R.A. and drawing a random decl. from the list of all ANTARES events. These coordinates serve as a sky position for the coincidence. The PSFs for $n_γ$ LAT photons and $n_ν$ neutrinos are then centered on this point, and placed randomly according to their respective PSFs. All photons are chosen to have the same inclination angle, which is drawn from the full set of inclination angles within the Fermi data set. A conversion type for each photon is similarly drawn from the Fermi data set. Photon energies are drawn from a power law with a photon index $Γ = 2$. Using the photon background maps, the number of unassociated photons expected to arrive within the temporal and spatial windows for that section of sky is calculated. From this Poisson probability, $n_b$ photons are randomly placed uniformly within the spatial window. Energy and conversion type for the background photons are chosen in the same manner as for the signal photons. All background photons are given the same inclination angle as the signal photons. Each particle is also given an arrival time randomly selected from a uniform distribution. Using this information, a $λ$ value is calculated following the methods of Section 3. Due to the iterative rejection of one or more low-significance γ-rays, events can end up with some of the injected photons excluded.

Because the varied physical models predicting $ν+γ$ coincidences have different characteristic timescales, we generate two sets of signal events for each of the three null distributions. One set draws the timestamps from a uniform distribution 100 s wide, while the other draws from a uniform distribution 1000 s wide.

To calculate the sensitivity of our analysis, we inject an increasing number of signal events $n_{inj}$ and plot the median resulting Anderson–Darling p-value (Scholz & Stephens 1987) against $n_{inj}/n_{obs}$ for the track and cascade data, as shown in Figure 3.

For the tracks, this provides an estimate of the threshold value of $n_{inj}$ that is needed to yield a statistically significant deviation from the null distribution (see columns $n_{inj,1%}$ and $n_{inj,0.1%}$ in Table 1). For the cascades, the size of each
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Anderson–Darling two-sample *p*-value vs. fraction of coincidences that result from signal events, *N*<sub>sig</sub>/N<sub>obs</sub>. Results from both signal populations are shown.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cumulative and residual histograms of the λ distributions for the track (left, $n_ν+γ = 2734$) and cascade (right, $n_ν+γ = 80$) data. The unscrambled data (green dashed line) and the null distribution (blue line) are shown for both tracks and cascades. Signal injections, generated using a 1000 s temporal window and yielding $p = 1\%$ (red line) and $p = 0.1\%$ (black line) are calculated for the track data only, as even 100% signal injection does not allow strong discrimination of signal and null distributions for the cascade data. Signal injection curves for the 100 s temporal window display as identical on this plot. Upper panels show cumulative histograms, while lower panels show residuals against the null distribution (plotted as null minus alternative). Anderson–Darling test *p*-values from comparison of the unscrambled and null distributions are $p = 39\%$ for the track sample and $p = 60\%$ for the cascade sample.</figcaption></figure>

<table><caption>**Table 2**
High-λ Events</caption><thead><tr><th>Date</th><th>Time (UTC)</th><th>MJD</th><th>Δt (s)</th><th>Position (J2000)</th><th>r<sub>1σ</sub></th><th>N<sub>ph</sub></th><th>λ</th><th>FAR (yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>2012 Nov 21</td><td>20:19:52</td><td>56252.8471</td><td>307</td><td>248°00, -7°70</td><td>2'</td><td>1</td><td>18.9</td><td>0.09</td></tr><tr><td>2014 Aug 5</td><td>11:13:33</td><td>56874.4677</td><td>750</td><td>279°68, -5°05</td><td>3'</td><td>2</td><td>18.8</td><td>0.09</td></tr></tbody></table>

**Note:** Date, Time, and MJD show the central time of the coincidence, while Δt measures the separation between the earliest and latest particles in the coincidence in seconds. Position gives the R.A. and decl. (in degrees) of the best-fit position, while $r_{1σ}$ gives the approximate $1σ$ error on the angular uncertainty in arcminutes (39% containment, assuming a Gaussian form). $N_{ph}$ is the number of photons in the coincidence. The false alarm rate (FAR) is calculated as the number of events of that λ or higher expected per year.

individual scramble is small enough that replacing 100% of the data set with signal events yields a *p*-value of 2.8% on average, making it very unlikely that this sample would yield a high-confidence demonstration of an underlying $ν+γ$ source population. At 90% confidence, our analysis is sensitive to >130 source-like $ν+γ$ coincidences in the 100 s track data, >145 in the 1000 s track data, and >60 in the 100 and 1000 s cascade data. Relevant statistics from these analyses are provided in Table 1.

In previous work, Turley et al. (2018) found that scrambled neutrinos coincident with LAT-detected GRBs, in particular GRB 090902B (Abdo et al. 2009), yielded λ values well above the $λ_C$ threshold. To quantify our analysis sensitivity to GRB + neutrino coincidences, we carried out a Monte Carlo simulation for each LAT-detected GRB<fnref n="57" /> that occurred within our data

<footnote n="57">⁵⁷ LAT GRB catalog: https://fermi.gsfc.nasa.gov/ssc/observations/types/grbs/lat_grbs/.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Average *Fermi* γ-ray background rates at the positions of track (upper panel) and cascade (lower panel) neutrinos. In each panel, the histogram shows the distribution obtained from 10,000 Monte Carlo scrambled data sets, while the red line marks the observed background rate for unscrambled data. Background rates are expressed in units of photons per square meter per square degree per 200 s. Observed average backgrounds are consistent with background for both data sets.</figcaption></figure>

collection period. Neutrinos were injected following our signal injection procedures, with the GRB position and trigger time as reference, and with a 1000 s box-window temporal distribution for neutrino arrival times. For each LAT GRB, we carried out 10,000 such neutrino signal injections and calculated the λ value for the resulting association in each instance.

The maximum λ generated through this search was λ = 3524.5, resulting from a 368-photon coincidence with GRB 130427A (Zhu et al. 2013). Of the 128 individual bursts in this simulation, 58 have median λ values from these neutrino injection trials of λ<sub>med</sub> > λ<sub>C</sub>, and a further five bursts have λ<sub>C</sub> > λ<sub>med</sub> > λ<sub>D</sub>.

4. Results

Applying our analysis to the two unscrambled neutrino data
sets yields the results summarized in Table 1. Figure 4 shows
the λ distributions for the unscrambled data for the track and
cascade data, along with the null distributions, and distributions
for signal injections (where possible) yielding *p*-values of 1%
and 0.1%, respectively.

All distributions are normalized to the number of coincidences in the unscrambled distribution. Note that due to the small size of the cascade coincidence sample, it is not possible to inject enough signal events into a random scramble to differentiate from other random scrambles at better than *p* = 2.8% (97.2% confidence).

Two coincidences above the λD threshold were observed in the track data. From Poisson statistics, two or more such coincidences would be observed 16.6% of the time given the 7.3 yr span of the data. Details of these two coincidences are presented in Table 2. No λ values above the λD threshold were observed in the cascade data. The subthreshold population search demonstrated that both unscrambled distributions were consistent with background, with test statistics of 39% for the tracks, and 60% for the cascades. Results from the track multiplet analysis are not shown as there were, on average, only 0.48 such coincidences per scramble, and none in the unscrambled analysis.

Turley et al. (2018) also tested for correlation between neutrino and *Fermi* LAT photon sky positions without any temporal correlation. Repeating this analysis using the ANTARES data, we first construct a single *Fermi* background map covering the full energy range. We then measure the background value at the location of every neutrino in the track and cascade data to compute an average photon background for each neutrino map. Carrying this out on the scrambled neutrino data sets yields an average background of (2.33 ± 0.06) × 10⁻² photons deg⁻² m⁻² per 200 s for the track data, and (2.16 ± 0.36) × 10⁻² photons deg⁻² m⁻² per 200 s for the cascade data. The observed backgrounds (in the same units) from the unscrambled data are 2.36 × 10⁻² (+0.44σ; *p* = 33%) for the track data, and 2.19 × 10⁻² (+0.09σ; *p* = 46%) for the cascade data. Both results are consistent with background (Figure 5).

The dispersion in the cascade background from scrambled
data sets is far larger than that for the tracks because of the
much-reduced sample size (180 cascade events compared with
7622 track events); however, the two average backgrounds are
consistent, as the mean of the track background is 0.47σ larger
than the mean of the cascade background, as measured using
the standard deviation of the cascade background distribution.
Recalling the IC 59 Northern (p = 28.1%), IC 59 Southern
(p = 4.7%), and IC 40 (p = 58.3%) results from Turley et al.
(2018), we can calculate a unified p-value of 19.7% from these
values using Fisher's method (Mosteller & Fisher 1948).

5. Conclusions

We have carried out a search for ν+γ transients using publicly
available *Fermi* LAT γ-ray data and ANTARES neutrino data.
Our analysis used archival data from both observatories over the
period 2008 August to 2015 December. As with previous work
(Turley et al. 2018), our analysis was designed to be capable of
identifying either individual high-significance ν+γ transients or a
population of individually subthreshold events, via statistical
comparison to uncorrelated (scrambled) data sets.

Our Monte Carlo simulations demonstrate a sensitivity to
single-neutrino events of sufficient γ-ray multiplicity, as demon-
strated by signal injection against multiple bright LAT-detected
γ-ray bursts. Signal injection against scrambled data sets
established our sensitivity to subthreshold populations of transient
ν+γ sources at the >7% level (>200 coincidences) for tracks;
however, due to the small sample size, we were not able to place
meaningful limits on a subthreshold ν+γ source population
within the cascades data. Our limit of >200 coincidences in the
full data set is equivalent to >27 LAT-associated cosmic
neutrinos per year in the ANTARES data. Because IceCube
estimates of the cosmic neutrino flux and spectrum lead us to
expect 6.8 cosmic ANTARES neutrinos per year (Section 1), our
limit is not physically constraining in this context.

Analysis of the observed (unscrambled) data reveals two
ν+γ coincidences above a nominal λD threshold (false alarm
rate FAR < 0.1 yr⁻¹; Table 2). Due to the 7.3 yr span of the
data, we anticipate observing two or more λ > λD coincidences
16.6% of the time (p = 16.6%). We observe no statistically
significant deviation of the observed λ distributions from their
associated null distributions, with observed *p*-values of
*p* = 39% and *p* = 60% for the track and cascade events,
respectively.

Independently, we performed the first test for correlation
between ANTARES neutrino positions and persistently bright
portions of the *Fermi* γ-ray sky. Our test found no significant
excess in either the tracks (p = 33%) or cascades (p = 46%). Combining these values with previous results (28.1% for IC 59 north, 4.7% for IC 59 south, 58.3% for IC 40; Turley et al. 2018) by Fisher's method yields a joint p-value of p = 19.7%.

While our results show no significant evidence of $\nu+\gamma$ coincidences, we look forward to the results of future searches using additional neutrino data. We also continue our work with Astrophysical Multimessenger Observatory Network (Smith et al. 2013; Cowen et al. 2016) partner facilities and the Gamma-ray Coordinates Network (Barthelmy et al. 1998) to generate low-latency $\nu+\gamma$ alerts from Fermi LAT $\gamma$-ray and high-energy neutrino data. Once these alerts are deployed, they will be distributed in real time to AMON follow-up partners.

The authors thank David Thompson for helpful discussions. We gratefully acknowledge support from Penn State's Office of the Senior Vice President for Research, the Eberly College of Science, and the Penn State Institute for Gravitation and the Cosmos. This work was supported in part by the National Science Foundation under Grant No. PHY-1708146. K.M. is supported by the Alfred P. Sloan Foundation and by the National Science Foundation under Grant No. PHY-1620777. 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), IdEx program and UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-INDEX-0005-02), Labex OCEVU (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-INDEX-0001-02), 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 Economía y Competitividad (MINECO): Plan Estatal de Investigación (refs. FPA2015-65150-C3-1-P, -2-P and -3-P, (MINECO/FEDER)), Severo Ochoa Centre of Excellence and Red Consolider MultiDark (MINECO), and Prometeo and Grisolía programs (Generalitat Valenciana), 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.

**Software:** Astropy (The Astropy Collaboration 2018), Matplotlib (Hunter 2007), HEASoft (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc), 2014), HEALPix (Górski et al. 2005), SciPy (Jones et al. 2001–2019)

## ORCID iDs

J. J. DeLaunay https://orcid.org/0000-0001-5229-1995

D. B. Fox https://orcid.org/0000-0002-3714-672X

A. Keivani https://orcid.org/0000-0001-7197-2788

K. Murase https://orcid.org/0000-0002-5358-5642

C. F. Turley https://orcid.org/0000-0002-9689-8075

A. Albert https://orcid.org/0000-0003-0197-5646

B. Baret https://orcid.org/0000-0001-6064-3858

S. Basa https://orcid.org/0000-0002-4291-333X

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

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

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

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

U. Katz https://orcid.org/0000-0002-7063-4418

R. Le Breton https://orcid.org/0000-0001-8522-4983

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

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

M. Mele https://orcid.org/0000-0002-9165-4231

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

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

## References

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

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017b, APh, 92, 30

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJL, 706, L138

Ackermann, M., Ajello, M., Allafort, A., et al. 2013, ApJ, 765, 54

Adrián-Martínez, S., Ageron, M., Albert, A., et al. 2016, JCAP, 2016, 062

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

Ageron, M., Aguilar, J. A., Al Samarai, I., et al. 2012, APh, 35, 530

Albert, A., André, M., Anghinolfi, M., et al. 2017a, EPJC, 77, 419

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

Albert, A., André, M., Anghinolfi, M., et al. 2017c, PhRVD, 96, 062001

Albert, A., André, M., Anghinolfi, M., et al. 2017d, ApJL, 850, L35

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

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

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

Albert, A., André, M., Anghinolfi, M., et al. 2019a, MNRAS, 482, 184

Albert, A., André, M., Anghinolfi, M., et al. 2019b, ApJ, 870, 134

ANTARES Collaboration 2015, JCAP, 12, 014

ANTARES Collaboration, Aguilar, A., Albert, J. A., et al. 2005, APh, 23, 131

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Barthelmy, S. D., Butterworth, P., Cline, T. L., et al. 1998, in AIP Conf. Ser. 428, Fourth Huntsville Gamma-ray Burst Symp., ed. C. A. Meegan, T. M. Koshut, & R. D. Preece (Melville, NY: AIP), 99

Cowen, D. F., Keivani, A., & Tešić, G. 2016, EPJ Web Conf., 116, 10001

Croft, S., Kaplan, D. L., Tingay, S. J., et al. 2016, ApJL, 820, L24

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hunter, J. D. 2007, CSE, 9, 90

IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2017, arXiv:1710.01191

IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018a, Sci, 361, 147

IceCube Collaboration, Fermi-LAT, MAGIC et al. 2018b, Sci, 361, eaat1378

Jones, E., Oliphant, T., Peterson, P., et al. 2001–2019, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/

Keivani, A., Fox, D. B., Tešić, G., Cowen, D. F., & Fixelle, J. 2015, arXiv:1508.01315

King, I. 1962, AJ, 67, 471

Kopper, C., Blaufuss, E. & IceCube Collaboration 2017, GCN Circ. 2196, https://gcn.gsfc.nasa.gov/gcn3/21916.gcn3

Mosteller, F., & Fisher, R. A. 1948, Am. Stat., 2, 30

Murase, K. 2015, AIP Conf. Proc., 1666, 040006

Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Scholz, F. W., & Stephens, M. A. 1987, J. Am. Stat. Assoc., 82, 918

Smith, M. W. E., Fox, D. B., Cowen, D. F., et al. 2013, APh, 45, 56

Tanaka, Y. T., Buson, S., & Kocevski, D. 2017, ATel, 10791, http://www.astronomerstelegram.org/?read=10791

The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Turley, C. F., Fox, D. B., Keivani, A., et al. 2018, ApJ, 863, 64

Turley, C. F., Fox, D. B., Murase, K., et al. 2016, ApJ, 833, 117

Zhu, S., Racusin, J., Kocevski, D., et al. 2013, GCN Circ. 14471, https://gcn.gsfc.nasa.gov/other/130427A.gcn3