Search for Gamma-Ray and Neutrino Coincidences Using HAWC and ANTARES Data

In the quest for high-energy neutrino sources, the Astrophysical Multimessenger Observatory Network (AMON) has implemented a new search by combining data from the High Altitude Water Cherenkov (HAWC) observatory and the Astronomy with a Neutrino Telescope and Abyss environmental RESearch (ANTARES) neutrino telescope. Using the same analysis strategy as in a previous detector combination of HAWC and IceCube data, we perform a search for coincidences in HAWC and ANTARES events that are below the threshold for sending public alerts in each individual detector. Data were collected between July 2015 and February 2020 with a livetime of 4.39 years. Over this time period, 3 coincident events with an estimated false-alarm rate of $<1$ coincidence per year were found. This number is consistent with background expectations.


ABSTRACT
In the quest for high-energy neutrino sources, the Astrophysical Multimessenger Observatory Network (AMON) has implemented a new search by combining data from the High Altitude Water Cherenkov (HAWC) observatory and the Astronomy with a Neutrino Telescope and Abyss environmental RESearch (ANTARES) neutrino telescope. Using the same analysis strategy as in a previous detector combination of HAWC and IceCube data, we perform a search for coincidences in HAWC and ANTARES events that are below the threshold for sending public alerts in each individual detector. Data were collected between July 2015 and February 2020 with a livetime of 4.39 years. Over this time period, 3 coincident events with an estimated false-alarm rate of < 1 coincidence per year were found. This number is consistent with background expectations.

INTRODUCTION
The Astrophysical Multimessenger Observatory Network (AMON Ayala Solares et al. 2020) is a virtual hub that integrates heterogeneous data from different astrophysical observatories with the main objective of enabling multimessenger astrophysics. Observatories that become members of AMON can act as trigger observatories or as follow-up observatories. Triggering observatories have high-duty cycles and a large field of view. Follow-up observatories have better angular resolution and sensitivity. AMON has developed coincidence analyses between high-energy gamma-ray and high-energy neutrino data. AMON mainly, but not necessarily, receives and uses data that are below the astrophysical event selection threshold (called subthreshold) for the individual observatories. In these data, possible signal events of astrophysical origin can be present and due to the limited sensitivity of a given detector (e.g., HAWC or ANTARES), cannot be distinguished from background events. Using a statistical analysis, AMON looks for temporal and/or spatial coincidences between events collected by different observatories with the purpose of recovering the signal events that are buried in the background.
The AMON analyses using gamma-ray and neutrino data include the coincidence studies between IceCube and Fermi -LAT (Turley et al. 2018); ANTARES and Fermi -LAT(Ayala Solares et al. 2019); and HAWC and IceCube (Ayala Solares et al. 2021). The last two analyses make use of the Neutrino-Electromagnetic (NuEM) AMON channel. This channel generates alerts in real time after receiving data from the respective observatories and performing a calculation to rank the coincidences (see Section §3). AMON servers are now located at the Amazon Web Services (AWS), having a high up-time (>99.99%). The NuEM alerts are sent as notices and circulars to the Gamma-ray Coordinates Network (GCN; Barthelmy 1990). Recently, AMON also started to send alerts to the Scalable Cyberinfrastructure to support Multi-Messenger Astrophysics (SCiMMA; Chang et al. 2019;Brady et al. 2019), a new hub for multimessenger astrophysics designed for private and public communication.
The NuEM channel searches for sources that emit secondary neutrinos and gamma rays. These neutrinos and gamma rays are produced in hadronic interactions, such as inelastic collisions of cosmic rays with matter or with radiation fields. These hadronic interactions produce neutral and charged pions, which then decay into the aforementioned particles. These interactions can occur in a wide variety of sources such as blazar flares, tidal disruption events, long gamma-ray bursts, short gamma-ray bursts, supernovae, and compact binary mergers (for a review of multimessenger sources, see Murase & Bartos 2019). In this work, we present a new analysis of this channel: the coincidence search between events collected by the HAWC gamma-ray observatory and the ANTARES neutrino telescope.
With ANTARES recently ceasing operations, this analysis helps us not only to look for possible sources in existing data, but also to prepare the necessary analysis tools and infrastructure for the KM3NeT neutrino telescope (Adrián-Martínez et al. 2016), the successor of ANTARES. The High Altitude Water Cherenkov (HAWC) observatory monitors the gamma-ray sky from its location in Puebla, Mexico, at an altitude of 4100 meters above sea level. Sitting between the volcanoes Sierra Negra and Pico de Orizaba, it has a large field of view that covers two-thirds of the sky daily. With a duty cycle above 95%, HAWC can monitor 2 sr of the sky continuously, which makes it ideal for observing transient events (Abeysekara et al. 2017). HAWC is a water Cherenkov detector array that characterizes the footprints of extensive air showers. Hadron-like showers and gamma-like showers can be classified by looking at how smooth and compact is the distribution of the charge measured by the photomultiplier tubes (PMTs) in the array. Hadron-like showers tend to have a discontinuous profile on the array due to the large number of muons in the shower, while gamma-ray showers present a smoother profile. By using the trigger time information of each PMT, reconstruction algorithms can find the direction of the primary particle with a 68% resolution of ∼ 0.2 • at energies above 10 TeV. HAWC is sensitive to gamma rays with energy from 300 GeV up to >100 TeV (Abeysekara et al. 2017;Albert et al. 2020).
The data that AMON receives from HAWC for this analysis includes the rising and setting time of the event position in the sky with respect to the detector -which defines the "HAWC transit time" of the event; a parameter (referred to as the "significance value" in the following) that estimates how much the event deviates from the expected hadron-like background and it is calculated after one transit; and the position in the sky of every event with their uncertainty. HAWC events are referred as "hotspots". The data used in this work were collected from July 2015 to February 2020.

ANTARES
The ANTARES neutrino telescope (Ageron et al. 2011) is located 40 km off-shore from the city of Toulon, France, in the Mediterranean Sea. It is a deep-sea Cherenkov neutrino detector. The detector consists of a three-dimensional array of 885 optical modules, each one with a 10 inch PMT, and distributed over 12 vertical strings anchored in the seafloor at a depth of about 2400 m. The detection of light from up-going charged particles is optimized with the PMTs facing 45 • downward. Since May 2007, the telescope has detected neutrino-induced muons that cause the emission of Cherenkov light in the detector, producing 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 . For the analysis presented in this manuscript, we use track-like events that are used in the point-source search analysis of ANTARES (Illuminati et al. 2019), which have a median angular resolution of 0.4 • for energies above 10 TeV. Since the ANTARES data is public, we are not using subthreshold data from ANTARES for this analysis.
The ANTARES data information consists of the following: the position and uncertainty of the individual observed event, the time of the event, and a p-value that quantifies the probability of the event to be a background event. For this study, we use the archival public data (2007-2017) that can be found in (Albert et al. 2018) as well as 3 more years of archival data (up to 2020) given by ANTARES through the AMON memorandum of understanding. Since we are using the ANTARES public data, it does not contain subthreshold events. We use the data that overlap with the time period of the HAWC data.

Computing the ranking statistics RS
The analysis method applied in this work is the same as the one developed for the HAWC and IceCube detectors in Ayala Solares et al. (2021), which is summarized below.
We assume to have a coincidence when the time of the ANTARES event falls between the rising and setting time of the HAWC event and the distance between the reconstructed directions of the events is smaller than 3.5 • (same as in Ayala Solares et al. 2021). After finding a coincidence, a test statistic is calculated to rank the coincidence. This is defined as which is based on the Fisher's method (Fisher 1938). The number of degrees of freedom of the test statistic is twice the number of p-values. The p λ value measures how much the events spatially overlap with each other. This value is obtained after optimizing a log-likelihood function defined as Gaussian on a sphere with x i and σ i being the measured position and positional uncertainty of the i-th event. B i is the background directional probability distribution from the corresponding detector at the position of the events. The sum is over all the N events that are part of the coincidence. The position of the coincidence, x coinc , is defined as the position where the log-likelihood is maximized, λ max . The p λ is obtained from the λ max distribution, which is the probability of seeing a λ max or higher.
The p HAWC value is related to the significance value of the HAWC event and quantifies the probability for the event to be from background.
The p Cluster is the probability of having n ν ANTARES events when one is already observed. It is defined as where ∆t is the HAWC transit time; f ν is the ANTARES background rate in a 3.5 • circle in the sky estimated as where f all is the measured background rate from the whole sky. Finally p ANTARES,i , is the fraction of ANTARES events that have a larger number of hits in the detector than the observed number of hits for the event. It is computed by using the normalized anti-cumulative distribution of the number of hits from the full ANTARES public data.
Since there can be n ν ANTARES events passing the selection criteria during a HAWC time window, the degrees of freedom of Eq. 1 vary. Therefore, we compute the p-value of the χ 2 6+2nν with 6 + 2n ν degrees of freedom. The ranking statistic (RS) is then simply defined as RS = − log 10 (p−value χ 2 6+2nν ).

Calculating the False Alarm Rate
The distribution of the RS is used to quantify the probability that the coincidences are fortuitous. It is also used to calculate the false alarm rate (FAR) of the coincidence (i.e. how rare the RS is). To build the distribution, we perform several simulations by scrambling the data sets a thousand times. The scrambling consists of randomizing the right ascension and the time values of the events. We then count how many events are above an RS value and divide by the simulated time (∼4600 years). Figure 1 shows the FAR as a function of the RS. A red line is shown which fits the data points, together with the 1σ uncertainty band.

Sensitivity and Discovery Potential
We obtain the sensitivity and discovery potential in the parameter space composed of the local rate density of transient sources vs the total neutrino isotropic energy as in Ayala Solares et al. (2021). We compare this to different source populations. We use the FIRESONG software package (Taboada et al. 2017) to obtain the number of transient sources during the same time as the archival data, along with the redshift, the neutrino flux normalization and the position in the sky of each source. The star formation rate assumed is from the core-collapse supernova rate obtained from the CANDELS and CLASH supernova surveys (Strolger et al. 2015). We assume a local rate density and a total neutrino isotropic equivalent energy as denoted along the x-and y-axes of Figure 2. The duration of each burst is fixed to 6 hours. For the neutrino energy spectrum, we assume a power law with a spectral index of -2.0 in the energy range between 10 TeV and 10 PeV. We use the model from Ahlers & Murase (2014) given as where d is the distance to the source; λ γγ is the interaction length of gamma rays with radiation backgrounds; K = 1 for photo-hadronic interactions and K = 2 for hadro-nuclear interactions; the sum is over the neutrino flavors. Using the neutrino flux normalization, and assuming photo-hadronic interactions, we can obtain the gamma-ray flux normalization from Eq. 5 1 . After obtaining the gamma-ray flux normalization, we inject the sources in HAWC simulated data. Using the simulated redshift information, we apply the attenuation of gamma rays from the extragalactic background light using the model from Domínguez et al. (2011). After running the HAWC analysis, if the observed hotspot has a significance larger than 2.75σ, we proceed to inject the neutrinos using Monte Carlo signal data from the ANTARES simulation as well as background events from the ANTARES scrambled data sets. Then we proceed to calculate the RS as explained in §3.1. We simulate transient sources for a period with the same livetime as the archival data being used in this analysis.
To be able to estimate the sensitivity and discovery potential, we use the number of coincidences above the 1 per year threshold as a statistic. For the livetime of the analysis, we expect to observe at least ∼ 4 random coincidences. Using random samples from the RS distribution, we find that the distribution of the number of coincidences behaves as a Poisson distribution with a mean of λ bk = 4.39, since that is the livetime of the analysis. We now need to find limits on the total signal and background rate , λ bk + λ s . In order to obtain the sensitivity we need a total rate that will produce a Poisson distribution with 90% of its population above the median of the background Poisson distribution. For the discovery potential, we need a total rate that will produce a Poisson distribution with 50% of its population above the threshold of the p-value = 2.87 × 10 −7 (5σ) of the background distribution. The mean signal λ s values for the sensitivity and discovery potential are 3.6 and 14.3.
We perform the simulation 100 times for each pair of local rate density of transient sources and total isotropic equivalent energies to gather enough statistics to build the Poisson distributions. The value of the rate densities and isotropic energies that give the desired λ s values are shown in Figure 2. The figure also shows several populations of transient sources, with a range of local rate densities and isotropic energies obtained from Murase & Bartos (2019). We see that long gamma-ray bursts are the only sources from which we may expect some detectable coincidences. Figure 2. Sensitivity, discovery potential (5σ) and 90 % upper limit for the archival data (analysis livetime of 4.39 years) in terms of total isotropic equivalent neutrino energy as a function of the local rate density. We assume a burst time of 6 hours and the neutrino spectrum to be a power law with index −2.0. Luminosity and rate-density ranges of the different sources can be found in Murase & Bartos (2019). For comparison, we show the sensitivity and discovery potential of the HAWC-IceCube analysis from Ayala Solares et al. (2021). Both the ANTARES effective area and the overlap region between HAWC and ANTARES are smaller compared to that in the HAWC-IceCube analysis. We can see that long gamma-ray bursts are potential candidates for a possible coincidence detection.

ARCHIVAL RESULTS
The ranking distribution for the archival data is shown in Figure 3, along with the simulated distribution from the scrambled data set used in Section 3.2 to obtain the FAR. The power of the combined data analysis can be seen in Figure 4, which shows the FAR for the HAWC-ANTARES coincidences versus the FAR for the HAWC events alone. As expected, the FAR for HAWC only events is reduced by 4 orders of magnitude. The low FAR of coincidences makes it useful for follow-up searches in real time.
After performing the analysis on unscrambled data, we found 3 events that pass the 1 per year FAR threshold in this period. These three events are summarized in Table 1. Although these events are rare given their FAR, they are still consistent with background as shown by the post-trials p-value (last column in Table 1, calculated as p post−trials = 1 − exp(−FAR · ∆T ), where ∆T is the livetime of the analysis). Tables 2 and 3 have the information of the events that make the coincidences.
The sky maps of the coincidences are shown in Figure 5. Each sky map shows the position of the individual events along with their uncertainties, as well as the best position of the coincidence. Also shown are the sources of the 4FGL Catalog (Abdollahi et al. 2020) that appear in each of the regions.   Table 1. Summary information for the three coincidences with FAR < 1 yr −1 . Information for the HAWC and ANTARES events that make these coincidences are found in Tables 2 and 3. The coincidence positions and uncertainty circles are shown as red in the skymaps of Figure 5. Note. The uncertainty U C (50%) corresponds to the 50% containment region of the estimated position of the coincidence. RS is the ranking statistic as defined in Section 3.1. The p-value corresponds to the post-trial p-value. Note. The uncertainty U H (50%) corresponds to the 50% containment region of the HAWC hotspot. The assumed flux model is a power law with an index shown in the last column. The index is fixed during the fit. The flux measurement is the normalization of the power law at 1 TeV. Note. The uncertainty U A (50%) corresponds to the 50% containment region of the ANTARES position. ∆θ is the distance from the best-fit HAWC hotspot position to the measured ANTARES event position.
We searched for past activity around the coincidences by looking in the Fermi All-sky Variability Analysis (FAVA) online tool 2 . We did not find any past activity in the regions of the coincidences 2 and 3. For coincidence 1, the source J0144.6+2705, associated with TXS 0141+268, is located 2.0 • away from the position of the coincidence. FAVA reported a burst in 2018 from this source, which was found in the high-energy band (800 MeV−300 GeV).
A search in the SIMBAD Catalog (Wenger et al. 2000) revealed several sources inside the uncertainty regions of the positions of the coincidences. For coincidence 1, we found several quasars, along with a radio source and an X-ray source. All of the quasars have redshift measurements larger than 0.3, the farthest HAWC can observe before the gamma rays start to be severely attenuated by the extragalactic background light (Albert et al. 2021).
In coincidence 2, there were 115 sources inside the uncertainty region of the coincidence in the SIMBAD Catalog. Around 15 sources are stars, while the rest are galaxies.
For coincidence 3 in 2019, we found only 8 sources in the SIMBAD catalog: 3 molecular clouds, 2 stars, 2 radio sources and one X-ray source. No information about the distance for the radio sources or X-ray source were available.
These coincidences are examples where, if the analysis had been running in real time, a follow-up observation in another wavelength could have pinpointed any source that is flaring in the region.

Upper Limit
After observing 3 coincidences in 4.39 years of data with a FAR of less than 1 per year, we calculate the 90% confidence level upper limit for the parameter space presented in Figure 2. Using Equation (9.54) from Cowan (2002), we obtain λ signal = 3.85. We apply the procedure of Section 3.3 to find the upper limit on the total isotropic equivalent neutrino energy as a function of the local source rate. (c) Coincidence 3. Figure 5. Sky maps in celestial coordinates of the HAWC-ANTARES coincidences with FAR values below 1 coincidence per year found in the archival data. The positions of the individual events are marked with the dots. The best-fit combined positions xcoinc, found after optimizing λ(x) (Eq. 2), are marked with a cross. Circles represent the 50% containment regions.

CONCLUSION
Archival data that span between 2015 and 2020 was analyzed to search for multimessenger sources through a coincidence analysis between subthreshold data of the HAWC observatory and public data from the ANTARES neutrino telescope. In this time period, three coincidences were found with a FAR of less than one coincidence per year. Although these coincidences are consistent with background expectations, they are still useful for follow-up observations, since the FAR can be improved by several orders of magnitude, compared to when the events are coming from the individual detectors. It is possible that a flare in the region could have been observed by a follow-up telescope, hinting at the presence of a multimessenger source, if the analysis had been running in real time.
Furthermore, based on the sensitivity and discovery potential studies, we found that long gamma-ray bursts are potential candidates for a possible coincidence detection with this analysis. This work is also a proof of principle analysis for future neutrino observatories. In this sense, with the end of operations of ANTARES, we expect that this analysis will be implemented for KM3NeT with an already exceeding ANTARES effective volume. Based on the information in Adrián-Martínez et al. (2016), a back-of-the-envelope estimate suggests an improvement of the sensitivity and discovery potential of more than one order of magnitude, assuming the same livetime. We look forward to implementing this new stream within AMON.