Multimessenger Gamma-Ray and Neutrino Coincidence Alerts Using HAWC and IceCube Subthreshold Data

The coincident detection of gravitational waves and electromagnetic radiation (Abbott et al. 2017), as well as the evidence found for a neutrino coincident with a gamma-ray flare from the blazar TXS 0506+056 (Aartsen et al. 2018), have shown the potential of multimessenger astrophysics. The ability to combine data from different observatories in real time or near-real time is driving this new era in astrophysics. The Astrophysical Multimessenger Observatory Network (AMON) has been created to facilitate the interaction of different observatories, create a framework for analyses with distinct data sets across multiple experiments, and notify the astrophysical community of any interesting events worthy of follow-up (Smith et al. 2013; Ayala Solares et al. 2020).⁸⁸

AMON focuses on using data that are below the discovery threshold of individual observatories. These events by themselves are heavily background-dominated, which complicates a search for astrophysical sources. By statistically combining the temporal and/or spatial information of these subthreshold events provided by different detectors, AMON aims to recover the signal events that are hidden among the background of each single observatory. Two multimessenger analyses were previously developed combining gamma-ray data from Fermi-LAT with neutrino data: one using IceCube data (Turley et al. 2018) and the other using ANTARES data (Ayala Solares et al. 2019).⁸⁹ The Fermi-LAT and ANTARES coincidence search started running in real time in 2019 April and has issued two alerts to date (see GCN circulars Turley 2020a, 2020b).

In this work, we focus on a new coincidence analysis combining information from the High Altitude Water Cerenkov (HAWC) Gamma-Ray Observatory (Abeysekara et al. 2017) and the IceCube Neutrino Observatory (Aartsen et al. 2017a) using the AMON infrastructure. This new multimessenger channel has been operational as a real-time coincidence search since 2019 December.

The purpose of this analysis is to search for hadronic accelerators that produce both gamma rays and neutrinos as secondary particles, with an emphasis on transient events. The accelerated cosmic rays can interact with target material surrounding the environment of the sources or with radiation fields. These interactions produce charged and neutral pions. Charged pions predominantly decay via π⁺ → μ⁺ + ν_μ, followed by the decay of the muon as ${\mu }^{+}\to {e}^{+}+{\nu }_{e}+{\bar{\nu }}_{\mu }$ (and charge conjugate). Neutral pions decay into two gamma-ray photons, π⁰ → γ + γ. The ratio between charged and neutral pions depends on the type of interaction of the cosmic rays with the targets. If the interaction occurs with electromagnetic radiation, the interaction will be photohadronic, which produces charged and neutral pions with probabilities of one-third and two-thirds, after considering both resonant and nonresonant pion productions. If the pions originate from interactions of cosmic rays with matter, the probability of producing charged and neutral pions is one-third for each type of pion (Biehl et al. 2019). A useful relation between the fluxes of gamma rays (F_γ) and neutrinos (${F}_{{\nu }_{\alpha }}$) is expressed as

$$\begin{eqnarray}&&{E}_{\gamma }{F}_{\gamma }({E}_{\gamma })\approx {e}^{-\tfrac{d}{{\lambda }_{\gamma \gamma }}}\displaystyle \frac{2}{3K}\displaystyle \sum _{{\nu }_{\alpha }}{E}_{\nu }{F}_{{\nu }_{\alpha }}({E}_{\nu }),\end{eqnarray} \tag{ 1 }$$

where E_γ ≈ 2E_ν are the gamma-ray and neutrino energies; α corresponds to the neutrino flavor; K is the ratio of charged to neutral pions, with K = 1 for photohadronic interactions and K = 2 for hadronuclear interactions; d is the distance to the source; and λ_γγ accounts for the attenuation of gamma rays due to their interaction with the extragalactic background light (EBL; see Murase et al. 2013 and Ahlers & Murase 2014).

In this paper, we present the algorithm and analysis to search for possible sources of gamma rays and neutrinos by looking at HAWC's and IceCube's subthreshold data. In Section 2, we briefly describe the detectors and their data. In Section 3, we present the statistical method and provide the false-alarm rate (FAR), sensitivities, and discovery potentials.

In Section 4, we present the results obtained using 3 yr of archival data, including upper limits for the same period of time for the total isotropic equivalent energy and source rate density parameter space. Finally, we conclude and discuss the implementation of the analysis in real time using the AMON framework.

The HAWC and IceCube are two detectors that focus on high-energy astrophysics, searching for sources that accelerate cosmic rays. Both detectors use the Cerenkov technique, where photomultipliers are used to detect the Cerenkov light produced by the passage of secondary charged relativistic particles—from gamma-ray, neutrino, and cosmic-ray showers—through a medium. The HAWC uses water as the medium, while IceCube uses the Antarctic ice.

Due to the attenuation of gamma rays on the extragalactic background photons, the signal from a source might not be significantly detected above background in the HAWC data. However, if IceCube neutrino events are found in spatiotemporal coincidence with a subthreshold HAWC hot spot, this might become an interesting coincidence that could be followed up by other observatories. In addition, the uncertainty region of HAWC events is generally smaller compared to IceCube events, which can give a better localization of a potential joint source.

2.1. High-energy Gamma Rays from HAWC

The HAWC observatory is a high-energy gamma-ray detector located in central Mexico. The complete detector has been in operation since 2015 March. The HAWC has a large field of view, covering two-thirds of the sky every day with a high-duty cycle in the decl. range from −26° to 64°. It is mainly sensitive to gamma rays in the energy range between 300 GeV and 100 TeV. It has an angular resolution of 02–10 (68% containment) that depends on the energy of the event, its zenith angle, and the size of the shower footprint measured by HAWC (Abeysekara et al. 2017).

We select locations of excess exceeding a given significance threshold—called "hot spots"—from the HAWC data to be used as inputs to the combined search. Hot spots are defined as locations in the sky with a cluster of events above the estimated cosmic-ray background level and measured by the significance (excess above the background). They are identified during one full transit of that sky location above the detector. The main hot-spot parameters AMON receives are the position coordinates and their uncertainty; the significance value, with a minimum of 2.75σ (threshold defined by HAWC); and the start and stop times of the transit. The duration of the transits are decl.-dependent, as shown in Figure 1. Since we are searching for unknown sources or sources that cannot be significantly detected above the background, we mask the data from the following parts of the sky above HAWC: the Galactic plane (b < ∣3°∣), the Crab Nebula, Geminga, Monogem, Mkr 421, and Mkr 501. The current rate of these hot spots received by AMON is ∼800 day^–1.

Zoom In Zoom Out Reset image size

2.2. High-energy Neutrinos from IceCube

The coincidence analysis is applied to events satisfying two criteria. The first is a temporal selection requiring the neutrino events to arrive within the transit time of the HAWC hot spot. Second, we select neutrinos that are within a radius of 35 from the HAWC hot-spot localization.⁹¹ After the neutrino events have passed the selection criteria, we calculate a statistic to rank the coincident events. The rate of coincidences after passing the criteria is ∼100 day^–1. This ranking statistic is based on Fisher's method (Fisher 1938), where we combine all the information that we have from the events. It is defined as

$$\begin{eqnarray}&&{\chi }_{6+2{n}_{\nu }}^{2}=-2\mathrm{ln}[{p}_{{}_{\lambda }}{p}_{{}_{\mathrm{HAWC}}}{p}_{{}_{\mathrm{cluster}}}\displaystyle \prod _{i}^{{n}_{\nu }}{p}_{{}_{\mathrm{IC},{\rm{i}}}}],\end{eqnarray} \tag{ 2 }$$

where the number of degrees of freedom is 6 + 2n_ν (as described below). The quantity ${p}_{{}_{\lambda }}$ quantifies the overlap of the spatial uncertainties of the events. The value ${p}_{{}_{\mathrm{HAWC}}}$ is the probability of the HAWC event being compatible with a background fluctuation. Since we can expect more than one IceCube candidate event in the time window (i.e., the HAWC transit period), we can calculate the probability of background IceCube events occurring in that time window. Given that we have at least one event detected, the ${p}_{{}_{\mathrm{cluster}}}$ ⁹² is the probability of that one event being in the same time window with the observed number of IceCube events, n_ν, or more from background; if there is only one IceCube event, this value is equal to 1.0. The value ${p}_{{}_{\mathrm{IC},{\rm{i}}}}$ is the probability of measuring a similar or higher energy/BDT score for an IceCube event, assuming it is a background event (calculated using the energy/BDT score and zenith angle). The ${p}_{{}_{\lambda }}$ value is obtained by a maximum-likelihood method that measures how much the positions of the HAWC and IceCube events overlap. This is calculated as

$$\begin{eqnarray}&&\lambda ({\boldsymbol{x}})=\displaystyle \sum _{i=1}^{N}\mathrm{ln}\left(\displaystyle \frac{{S}_{i}({\boldsymbol{x}})}{{B}_{i}}\right),\end{eqnarray} \tag{ 3 }$$

where N is the HAWC hot spot plus the number of IceCube candidate events; S corresponds to a signal directional probability distribution function, which is assumed to be a Gaussian distribution on the sphere with a width given by the measured positional uncertainty from each detector, ${S}_{i}{({\boldsymbol{x}})=\exp [-({\boldsymbol{x}}-{{\boldsymbol{x}}}_{i})}^{2}/2{\sigma }_{i}^{2}]/(2\pi {\sigma }_{i}^{2})$; and B_i is the background directional probability distribution from each detector at the position of the events. The position of the coincidence, x_coinc, is defined as the position of the maximum-likelihood value, ${\lambda }_{\max }$, as shown in Figure 5. The uncertainty of x_coinc is calculated by the standard error ${\sigma }_{{{\boldsymbol{x}}}_{\mathrm{coinc}}}^{2}=1/{\sum }_{i}^{N}({\sigma }_{i}^{-2})$.

The λ_max values are used to make a distribution of the overlap of the coincidences. A higher λ_max value indicates a more significant overlap of the event uncertainties. This translates into a smaller p-value ${p}_{{}_{\lambda }}$.

Due to the fact that we can have more than one IceCube event passing the selection criteria, the degrees of freedom of Equation (2) vary. We therefore calculate a p-value of the χ² with 6 + 2n_ν degrees of freedom. The ranking statistic is then simply defined as $-{{\rm{l}}{\rm{o}}{\rm{g}}}_{10}$ (p-value).

3.1. Calibration of the FAR

We apply the above-described algorithm to 3 yr of scrambled data sets from both observatories. Scrambling consists of randomizing the R.A. and time values of the events many times in order to calibrate the FAR. The result of this process is shown in Figure 2. For a specific ranking statistic, we calculate the total number of coincidences above this ranking statistic value and then divide by the total amount of scrambled simulation time to get the rate. The linear fit in Figure 2 is used to estimate the FAR in real-time analyses.

Zoom In Zoom Out Reset image size

3.2. Sensitivity and Discovery Potential

To put the archival results into context, we look at a simulation for transient events that can produce both neutrinos and gamma rays. We quantify the sensitivity and discovery potential for the 1 coincidence yr^–1 threshold for a live time of 3 yr of data.

We use the FIRESONG software package (Taboada et al. 2018), which simulates neutrino sources for a given local rate density of transient gamma-ray and neutrino sources, total neutrino isotropic equivalent energies, and timescales. The outcome of the simulation is a list of simulated neutrino sources with decl., redshift, and neutrino flux normalization. This is based on a power-law energy spectrum with a spectral index of −2 for the flux, in the energy range between 10 TeV and 10 PeV,⁹³ and a time of the burst of 6 hr.⁹⁴ Using Equation (1), we can transform the normalization to a gamma-ray flux assuming photohadronic interactions. We then simulate the sources in HAWC, also adding EBL attenuation with the model from Domínguez et al. (2011), and in addition, we draw a Poisson random number of neutrinos with an expectation value given by the source flux and IceCube's background. Finally, we quantify the coincidence.

We calculate the sensitivity and discovery potential by running simulations for a given pair of rate density and total neutrino isotropic energy. We apply the coincidence algorithm and, after finding the signal coincidences, add them to a distribution with random coincidences. We keep the total number of coincidences the same as that of the 3 yr of data, so we remove the same number of random coincidences as injected sources. We apply this procedure several times in order to build a distribution of the number of coincidences that cross the 1 coincidence yr^–1 threshold, N(FAR ≤ 1). If no sources are injected, N(FAR ≤ 1) is a Poisson distribution with a rate of r_B = 3.0 (B stands for background) for the 3 yr of observations. For the sensitivity, we find the pair of parameters that will give us r_B + r_S = 6.0 (where S stands for signal). This corresponds to an N(FAR ≤ 1) distribution that crosses the median of the Poisson background distribution 90% of the time. For the 5σ discovery potential, we find the pair of parameters that will give a rate of r_B + r_S = 15.7, since this distribution will have 50% of its population with a p-value smaller than 2.87 × 10⁻⁷ with respect to the Poisson background distribution. We fit the distribution of N(FAR ≤ 1) to a Poisson function and find the best value for r_S. The pair of rate density and total neutrino isotropic energy that gives the corresponding r_S values for sensitivity or discovery potential is plotted in Figure 3. To put the sensitivity and discovery potential in context, we include diagonal lines that show the total neutrino isotropic energy as a function of the rate density that would be required to produce the total observed IceCube diffuse neutrino flux (assuming a power-law spectrum with an index of −2.5). This assumes either no evolution or star formation evolution following the Madau–Dickinson model (Madau & Dickinson 2014); it also assumes a standard candle luminosity function. Based on Aartsen et al. (2018), we marked a region in Figure 3 showing the estimated released neutrino energy of the IceCube event 170922A related to TXS 0506+056.

Zoom In Zoom Out Reset image size

4.1. Archival Data

We analyzed data collected from 2015 June to 2018 August. Figure 4 shows the distribution of the ranking statistic value of the unblinded data compared to the expected distribution of random coincidences (i.e., the scrambled data sets mentioned in Section 3).

Zoom In Zoom Out Reset image size

Since we are interested in searching for rare coincidences, we look for coincidences with an FAR of less than 1 coincidence yr^–1, which corresponds to a ranking statistic value of 7.31. We found two coincidences, one in 2016 and one in 2018, with ranking statistics of 7.34 (1 coincidence yr^–1) and 9.43 (1 coincidence in 38.5 yr), respectively. These coincidences are not significant with respect to the background distribution. Using p-value = $1-\exp (-t\cdot {\rm{FAR}})$, with t = 3 yr, the p-values are 0.95 and 0.075, respectively. The sky maps of the two coincident events with the highest ranking statistic values are shown in Figure 5. Table 1 contains the summary information on them. Information on the individual events that form each coincidence can be found in Tables 2 and 3.

Zoom In Zoom Out Reset image size

Table 1.  Summary Information on the Two Coincidences with FAR < 1

Decl. (deg)	R.A. (deg)	Uncertainty (50% Containment)(deg)	Ranking Statistic	FAR (yr–1)	p-value
2.96	4.93	0.16	7.3	0.99	0.95
2.27	173.99	0.53	9.4	0.026	0.075

Download table as:  ASCIITypeset image

Table 2.  Information on the Two HAWC Hot Spots that Correspond to Each of the Coincidences with an FAR < 1 yr^–1 in the 3 yr Data Set

Decl.	R.A.	Uncertainty	Initial Time	Final Time	Significance	Flux Upper Limit
(deg)	(deg)	(deg)	(UT)	(UT)	σ	(TeV−1 cm−2 s−1)
2.91	4.96	0.17	2016-12-11 22:11:47	2016-12-12 04:38:41	3.71	3.9e-11
2.38	173.4	0.74	2018-04-12 01:31:21	2018-04-12 07:54:51	2.77	8.3e-11

Note. Flux upper limits are based on an E⁻² energy spectrum.

Download table as:  ASCIITypeset image

Table 3.  IceCube Neutrino Information for Each of the Coincidences

Decl.	R.A.	Uncertainty	Time	Background p-value	Δθ
(deg)	(deg)	(deg)	(UT)	${p}_{{}_{{IC}}}$	(deg)
3.04	6.86	1.31	2016-12-11 23:20:25	0.944	1.90
2.66	4.35	0.71	2016-12-12 00:24:48	0.055	0.65
5.18	3.00	1.08	2016-12-12 01:37:28	0.391	2.99
5.71	6.92	2.13	2016-12-12 03:22:12	0.993	3.42
0.30	172.77	1.67	2018-04-12 01:57:33	0.222	2.12
4.45	174.88	1.61	2018-04-12 03:53:08	0.860	2.51
1.75	175.88	1.48	2018-04-12 04:36:11	0.001	2.50
2.05	174.42	1.42	2018-04-12 05:19:36	0.005	1.02

Note. The uncertainty corresponds to the 50% containment. Here Δθ is the distance from the best-fit HAWC hot-spot position to the measured neutrino position.

Download table as:  ASCIITypeset image

We looked at the SIMBAD catalog (Wenger et al. 2000) for sources that appear near the coincidences⁹⁵ and the Fermi All-sky Variability Analysis (FAVA) online tool⁹⁶ for any evidence of past flares in the region based on the light curves provided by FAVA.

For the coincidence of 2016 with an FAR of 0.99 yr^–1, there is a radio galaxy in the nearby region, PKS 0017+026, also known as TXS 0017+026 (Dunlop et al. 1989). This source is 004 away from the best-fit position of the coincidence. Unfortunately, no distance information is available to estimate the gamma-ray attenuation. Other sources that appear nearby are quasars, but in general, these sources are too distant (redshift above 0.3), resulting in strong gamma-ray attenuation. With the FAVA tool, the source from the 3FGL catalog, J0020.9+0323, was found 052 away from the best-fit coincidence position, which is outside the 50% containment region. The 3FGL catalog mentions that this is an unassociated source (Acero et al. 2015).

For the coincidence of 2018 with an FAR of 0.026 yr^–1, several sources appear in the SIMBAD catalog. There are nine radio galaxies within 074 of the best-fit location of the coincidence from the NRAO VLA Sky Survey Catalog. From these, only NVSS J113719+022200 had some information about its distance (redshift of 0.19). We did not find nearby sources in the FAVA monitoring tool for this coincidence.

Both coincidences found with this analysis are therefore consistent with background expectations. Follow-up observations in the optical and X-ray could be helpful to discern if any of these sources are related to the coincident events.

4.2. Upper Limit

The real-time implementation of the analysis started on 2019 November 20. As specified in Ayala Solares et al. (2020), we use the amonpy software for the real-time implementation of the analysis. A major difference is that the system is now running on Amazon Web Services servers, which will further improve AMON's uptime. We set a threshold for public alerts at an FAR < 4 coincidences yr^–1. This threshold is set so that there is a reasonable number of statistically interesting coincidences that can be followed up during a year. Alerts are sent immediately to AMON members, and a GCN notice is generated. A GCN circular is also written to inform the rest of the astrophysical community. The first public alert of the system was sent out on 2020 February 2. It had an FAR of 1.39 yr^–1. The reported position is (R.A., decl.) = 2003, 1271, with a 50% radius of 017 (see GCN circular 26963; Ayala Solares 2020). The MASTER Global Robotic Net and the ANTARES observatory performed follow-up observations of the coincidence, but no transient event was observed (see GCN circulars 26973 and 26976; Kouchner 2020; Lipunov 2020).

The largest latency of the analysis comes from the HAWC analysis of the hot spots, since the transit needs to complete before sending that information to AMON. Based on Figure 1, the hot-spot duration can last from less than 1 hr to a bit more than 6 hr. The latency, once the data are in the AMON server, is less than 1 minute to perform the analysis and send the alert to the public.

We developed a method to search for coincidences of subthreshold data from the HAWC and IceCube observatories. Using coincidences of subthreshold data allows us to recover signal events that cannot be differentiated from the background in each individual detector. The method was tested on archival data taken between 2015 and 2018. We found two coincidences in the archival analysis that crossed the FAR threshold of 1 yr^–1, consistent with the background expectations of three coincidences in 3 yr. Although a few sources were found near the best coincidence positions, these results are still consistent with the expectation from random coincidences. The real-time analysis has produced one alert so far, with an FAR of 1.39 yr^–1. It was sent out to the community. We encourage other observatories to perform follow-up observations of these results and the real-time alerts in the future.

AMON: This research or portions of it were conducted with Advanced CyberInfrastructure computational resources provided by the Institute for Computational and Data Sciences at the Pennsylvania State University (https://www.icds.psu.edu/). This material is based upon work supported by the National Science Foundation under grants PHY-1708146 and PHY-1806854 and the Institute for Gravitation and the Cosmos of the Pennsylvania State University. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. F.K. was supported as an Eberly Research Fellow by the Eberly College of Science at the Pennsylvania State University.

HAWC: We acknowledge support from the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants IG101320, IN111315, IN111716-3, IN111419, IA102019, IN112218; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society—Newton Advanced Fellowship 180385; Generalitat Valenciana, grant CIDEGENT/2018/034; and Chulalongkorn University's CUniverse (CUAASC) grant. Thanks to Scott Delay, Luciano Díaz, and Eduardo Murrieta for technical support.

IceCube: We gratefully acknowledge the following support: USA—US National Science Foundation Office of Polar Programs, US National Science Foundation Physics Division, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin–Madison, Open Science Grid (OSG), Extreme Science and Engineering Discovery Environment (XSEDE), US Department of Energy National Energy Research Scientific Computing Center, Particle Astrophysics Research Computing Center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle Physics Computational Facility at Marquette University; Belgium—Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programs, and Belgian Federal Science Policy Office (Belspo); Germany—Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing Cluster of the RWTH Aachen; Sweden—Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia—Australian Research Council; Canada—Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark—Villum Fonden, Danish National Research Foundation (DNRF), Carlsberg Foundation; New Zealand—Marsden Fund; Japan—Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea—National Research Foundation of Korea (NRF); Switzerland—Swiss National Science Foundation (SNSF); United Kingdom—Department of Physics, University of Oxford.

Facilities: HAWC - , IceCube - , AMON. -

Software: astropy (Price-Whelan et al. 2018), FIRESONG (Taboada et al. 2018), numpy (Van der Walt et al. 2011), scipy (Virtanen et al. 2020) matplotlib (Hunter 2007), pandas (McKinney 2010), amonpy (Ayala Solares et al. 2020).