A Coincidence Search for Cosmic Neutrino and Gamma-Ray Emitting Sources Using IceCube and Fermi-LAT Public Data

Our analysis begins by filtering for coincidences between an individual neutrino event and all photons within 5° angular separation and ±100 s arrival time, as per Keivani et al. (2015). The angular acceptance cut corresponds approximately to the maximum 1σ radial uncertainty for Fermi-LAT photons satisfying our event selection. The temporal acceptance window is chosen to include ≈90% of classical GRBs (Paciesas et al. 1999). For each coincidence, a pseudo-log-likelihood test statistic, λ, is calculated as follows:

$\begin{eqnarray}&&\lambda =2\mathrm{ln}\displaystyle \frac{({P}_{\gamma 1}({\boldsymbol{x}}){P}_{\gamma 2}({\boldsymbol{x}})...{P}_{\gamma n}({\boldsymbol{x}}))n!({P}_{\nu }({\boldsymbol{x}}))}{{B}_{1}({\boldsymbol{x}},{E}_{1},{\theta }_{1}){B}_{2}({\boldsymbol{x}},{E}_{2},{\theta }_{2})...{B}_{n}({\boldsymbol{x}},{E}_{n},{\theta }_{n})},\end{eqnarray} \tag{ 1 }$

where n is the number of photons coincident with the neutrino, ${P}_{\gamma i}({\boldsymbol{x}})$ is the energy-dependent PSF of the LAT at the best-fit position, ${\boldsymbol{x}}$, and ${P}_{\nu }({\boldsymbol{x}}))$ is the IceCube PSF at the best-fit position. Both PSFs have units of probability per square degree. The ${B}_{i}({\boldsymbol{x}},{E}_{i},{\theta }_{i})$ are the LAT background terms in units of photons per square meter (approximating the Fermi effective area) per 200 s (our temporal window) per square degree for each ${\gamma }_{i}$, given its energy E_i and inclination angle ${\theta }_{i}$. In this metric, larger values of λ indicate a higher likelihood correlated multiplet. This pseudo-log-likelihood statistic is the natural extension of the Keivani et al. (2015) test statistic to multi-photon coincidences, via the Poisson likelihood of generating an n-fold coincidence from background; in the prior approach, each ν + γ coincidence was treated separately.

The best-fit position ${\boldsymbol{x}}$ is determined as the location of maximal PSF overlap. As the overlap of a double King function with a Gaussian function cannot be solved analytically, the best-fit position is found numerically. For single neutrino + multiple photon coincidences, the event photon multiplicity is determined by optimization: we compare the λ value at maximum multiplicity to that which would result if the photon with the lowest PSF density at the best-fit position were excluded (after recalculating the best-fit position and λ), and iteratively exclude photons until λ no longer increases.

We generate a set of 10,000 Monte Carlo scrambled versions of each of our three data sets in order to characterize their null distributions and define analysis thresholds, prior to performing any analysis of the unscrambled data sets. Our scrambling procedure begins by shuffling the full set of neutrino detections, associating each original neutrino ${\nu }_{i}$ with another randomly selected neutrino ${\nu }_{j}$. Each neutrino ${\nu }_{i}$ retains its original decl. and angular error and receives the original arrival time of neutrino ${\nu }_{j}$, with its new R.A. derived by adjusting the original R.A. for the difference in local sidereal time between the original and new arrival times, the same approach as in Turley et al. (2016). 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).

This analysis presents two discovery scenarios. First, since our test statistic λ is unbounded, the null distribution provides us with threshold values, which can be used to identify individually significant coincidences and estimate their false alarm rates. We define two such thresholds, ${\lambda }_{\times 10}$, the value exceeded (one or more times) in 1 of 10 scrambled data sets, and ${\lambda }_{\times 100}$, the value exceeded in 1 of 100 scrambled data sets. For analyses treating a year of observations, these two thresholds would correspond to events with false alarm rates of 1 decade⁻¹ and 1 century⁻¹, respectively.

Under this approach (and without accounting for the trials factor, see below), observation of a single event above ${\lambda }_{\times 100}$, or two events above ${\lambda }_{\times 10}$, would constitute evidence of joint ν + γ emitting sources, while observation of two events above ${\lambda }_{\times 100}$ or four events above ${\lambda }_{\times 10}$ would enable a discovery claim.

In the absence of any individually significant events, there remains the opportunity for discovery of a subthreshold population of ν + γ emitting sources. By design, true cosmic coincidences are biased to higher λ values (Figure 3), and a population containing a sufficient number of such signal events can be distinguished from the null distribution using the Anderson–Darling k-sample test. The k-sample test is used to establish mutual consistency among k observed data sets (k = 2 for a two-sample test), testing against the null hypothesis that they are drawn from a single underlying distribution.

Zoom In Zoom Out Reset image size

Given our choice to make two statistical tests on each of three predefined data sets (IC 40, IC 59-North, and IC 59-South), we apply an ${N}_{\mathrm{trials}}=6$ trials penalty to our unscrambled analyses (Section 4.1).

To estimate our sensitivity to cosmic ν + γ emitting source populations, we generate signal-like events and inject these into scrambled data sets, comparing the results to the null distribution.

Since we test for γ multiplicity, as part of this process we must adopt a procedure for determining whether each injection is a γ singlet, doublet, or ${n}_{\gamma }\gt 2$ multiplet. To determine the appropriate n_γ distribution, we assume a population of sources emitting one neutrino, with associated photon fluence distributed according to $N(S\geqslant {S}_{0})\propto {S}_{0}^{-3/2}$, where S₀ is a threshold photon fluence and $N(S\geqslant {S}_{0})$ is the number of events observed with fluences greater than or equal to this threshold; we note that an ${S}_{0}^{-3/2}$ dependence is expected for source populations of arbitrary luminosity function distributed in Euclidean space.

Adopting a minimum-considered fluence of ${S}_{\min }=0.001$ photons and inverting this relation, we generate the expectation value for the multiplicity of any event as $\langle {n}_{\gamma }\rangle ={S}_{\min }\,{u}^{-2/3}$, where u is a uniform random variable. We then generate the observed n_γ by drawing randomly from the Poisson distribution with expectation value $\langle {n}_{\gamma }\rangle$. Excluding zero-multiplicity events, we are left with the following n_γ distribution: 93% singlet, 4.8% doublet, 1.1% triplet, 0.5% quadruplet, 0.3% quintuplet, 0.2% sextuplet, and 0.1% septuplet (the highest multiplicity we allow). As an aside, we note that this is (approximately) the unique and universal distribution expected to arise in these cases, for extragalactic source populations extending to modest redshift ($z\lesssim 1$) with weak evolution, and thus distributed in near-Euclidean space.

A signal event of photon multiplicity n_γ is generated by centering the PSF for n_γ LAT photons and an IceCube neutrino at the origin. The neutrino localization uncertainty is drawn from the full set of IceCube neutrino uncertainties, while the inclination angles and conversion types of the photons are drawn from the full set of these distributions within the Fermi data set. Photon energies are drawn from a power law with a photon index ${\rm{\Gamma }}=-2$. The photons and the neutrino are placed randomly according to their respective PSFs. A random sky position is then chosen as the best-fit position for this coincidence, and a λ value is calculated following the methods of Section 3. Since the λ calculation involves maximizing λ by exclusion of outlying photons, many events end up with some of the injected photons excluded. Cumulative distributions for the null and signal-only distributions are shown in Figure 3.

To calculate the sensitivity of our analysis, we inject an increasing number of signal events ${n}_{\mathrm{inj}}$ into a scrambled (null) distribution and compare the signal-injected and null λ distributions using the Anderson–Darling k-sample test. We carry out 10,000 trials for each selected value of ${n}_{\mathrm{inj}}$ and plot the mean resulting p-value against ${n}_{\mathrm{inj}}$, for each of our data sets, in Figure 4. In this way we estimate the threshold value of ${n}_{\mathrm{inj}}$ that is required to yield a statistically significant deviation from the null distribution for each of the data sets (${n}_{\mathrm{inj},1 \% }$ and ${n}_{\mathrm{inj},0.1 \% }$ columns in Table 1).

Zoom In Zoom Out Reset image size

Table 1.  Coincidence Search Results

		Thresholds				Observed
Data Set	$\langle {n}_{\nu +\gamma }\rangle$	${\lambda }_{\times 10}$	${\lambda }_{\times 100}$	${n}_{\mathrm{inj},1 \% }$	${n}_{\mathrm{inj},0.1 \% }$	${n}_{\nu +\gamma }$	${\lambda }_{\max }$	${p}_{{\rm{A}}-{\rm{D}}}$
IC 40	1090 ± 30	23.9	27.2	150	210	1128	20.3	63%
IC 59-North	4970 ± 65	26.5	49.0	440	570	5046	17.8	16.8%
IC 59-South	7072 ± 76	26.8	31.5	565	740	7080	24.4	3.8%

Note. $\langle {n}_{\nu +\gamma }\rangle$ is the expected number of neutrinos observed in coincidence with one or more gamma rays, as derived from 10,000 Monte Carlo scrambled realizations of each data set. ${\lambda }_{\times 10}$ and ${\lambda }_{\times 100}$ are the thresholds above which a coincidence is only observed once per 10 or 100 scrambled data sets, respectively. ${n}_{\mathrm{inj},1 \% }$ and ${n}_{\mathrm{inj},0.1 \% }$ are the number of injected signal events required in simulations to give an Anderson–Darling test statistic of $p\lt 1 \%$ and $p\lt 0.1 \%$, respectively, by comparison to the null distributions for each data set. ${n}_{\nu +\gamma }$ is the number of neutrinos observed in coincidence with one or more gamma rays in the unscrambled data, ${\lambda }_{\max }$ is the maximum observed λ for each data set, and ${p}_{{\rm{A}}-{\rm{D}}}$ is the value of the Anderson–Darling test statistic from comparison of the observed λ distribution to its associated null distribution.

Download table as:  ASCIITypeset image

Carrying out the scrambled analysis on the IC 40 and IC 59 data produces the null λ distributions shown in Figure 3. Key statistics from these analyses are summarized in Table 1. Most of the simulated events with $\lambda \gt {\lambda }_{\times 100}$ in IC 59-North scrambled runs result from a scrambled neutrino landing in near coincidence with one of two GRBs detected by the LAT during our period of observation. GRB 090902B (Abdo et al. 2009) placed over 200 photons on the LAT, giving a maximum $\lambda =2560.2$ in a 218-photon coincidence. GRB 100414A (Takahashi et al. 2010) placed over 20 photons on the LAT and yields a maximum $\lambda =91.2$ in a 10-photon coincidence. Excluding all coincidences with either of these GRBs would yield a threshold of ${\lambda }_{\times 100}=35$ for the IC 59-North data, rather than the GRB-inclusive value of ${\lambda }_{\times 100}=49.0$.

Given the number of signal-like ν + γ required to yield a $p\lt 1 \%$ deviation in the Anderson–Darling k-sample test, we estimate our analysis would detect >150 source-like ν + γ coincidences from IC 40 (>13% of the total number of coincidences in IC 40), >440 from IC 59-North (>9%), and >565 from IC 59-South (>9%); see Table 1.

Applying our analysis to the three unscrambled neutrino data sets yields the results summarized in Table 1. Figures 5 and 6 show the λ distributions for the unscrambled data for IC 40, IC 59-North, and IC 59-South, along with the null distributions, and distributions for signal injections yielding p-values from the Anderson–Darling test of 1% and 0.1%, respectively. All distributions are normalized to the number of coincidences ${n}_{\nu +\gamma }$ observed in the unscrambled data. No λ values were detected above the ${\lambda }_{\times 10}$ threshold in any of the analyses. Notably, as seen in Figure 6, the IC 59-North and IC 59-South data show an excess of lower λ values by comparison to the null distributions, unlike the excess of higher λ values expected from a signal population.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given our six trials, a minimum observed single-trial $p=3.8 \%$ corresponds to a trials-corrected value of ${p}_{\mathrm{post}}=20.7 \%$, which is not significant. However, the similar scale and shape of the residual patterns for IC 59-North and IC 59-South lead us to seek out possible causes of these residual patterns. To illustrate this point, combining p-values from these two data sets (our most sensitive) by Fisher's method gives a joint $p=3.9 \%$ (single trial) as the probability of generating two such large deviations by random chance.

We note that we have not conceived of any way for systematic effects to generate the IC 59 residuals and low p-values, simply because any errors or simplifications in the analysis (which certainly exist) are replicated across all scrambled data sets. Rather, the only way to generate these effects (if they are not due to random statistical fluctuation) is via spatio-temporal correlation of neutrinos and gamma rays. Such correlation would imply either cosmic sources, or at a minimum, correlated emission (e.g., enhancement toward the Galactic plane or Supergalactic plane) and hence, require structure in the neutrino sky, which has not previously been observed.

We divide our further explorations into two approaches: first (Section 4.2), we further vet the IC 59 data against our original hypothesis, to test for any evidence that short-duration ($\delta t\lt 100$ s) ν + γ transients are really present in the data. Second (Section 4.3), we test for longer-duration ν + γ spatio-temporal correlations that might have an effect on our original analysis.

We vet the IC 59 data sets to evaluate whether the low p-values and systematic trends in the residual λ distributions that we observe could be due to ν + γ transient sources, as per our original hypothesis, but below the sensitivity of those analyses. We exclude the IC 40 data set from these tests as it is less sensitive to the presence of cosmic sources (see Figure 4 and Section 3.4).

Specifically, we check for systematic trends or anomalies in the spatial and temporal distributions of the neutrino-coincident photons that might account for the unexpected deviation to lower λ values. We test separately for deviations in the distributions of the photons' angular and temporal separations from their coincident neutrino, for IC 59-North and IC 59-South.

With regards to the angular separation distributions, we note that a systematic underestimation of IceCube neutrino localization uncertainties might cause suppressed λ values relative to simulations, due to the P_ν term in the pseudo-likelihood. In a similar vein, even if all neutrino and gamma-ray localization uncertainties are accurately characterized, a systematically softer spectrum for cosmic ν + γ sources (${\rm{\Gamma }}\lt -2$) would suppress λ values via the ${P}_{\gamma i}$ terms in the pseudo-likelihood, since higher-energy LAT photons are better localized.

We construct five-bin histograms of the angular separations of all neutrino-associated photons, with bin boundaries chosen to make the null distribution approximately flat (equal numbers of photons in each bin). We then calculate the ${\chi }^{2}$ statistic for the unscrambled data compared to the (flat) null distribution.

Figure 7 (note zero-suppressed y axis) presents our results. Observed IC 59-North angular separations are consistent with the (flat) null distribution, exhibiting ${\chi }^{2}=8.238$ for 4 degrees of freedom ($p=14.3 \%$). IC 59-South angular separations, by contrast, show a substantial deficit in the bin at $\delta \theta \approx 3^\circ$ (and modest excesses in the bins to either side), which results in ${\chi }^{2}=11.502$ for 4 degrees of freedom, giving $p=4.2 \%$. While this deviation is moderately surprising, the absence of any systematic trend to low or high separations suggests it is likely not responsible for the observed deviation in the λ distribution.

Zoom In Zoom Out Reset image size

Here we note that a systematic trend to small angular separations would suggest the presence of ν + γ sources as per our test hypothesis, while a systematic trend to large angular separations would suggest the presence of ν + γ sources with underestimated localization uncertainties or soft gamma-ray spectra. Neither such trend is observed.

We execute a similar analysis of the temporal separations of coincident photons. In contrast to the angular separations, which are incorporated into our pseudo-likelihood calculation, temporal separations are not considered (apart from the predefined acceptance window), so this analysis serves as an independent test of our original hypothesis. For purposes of trials correction of any subsequent statistics, we therefore add two trials, giving ${N}_{\mathrm{trials}}=8$.

Figure 8 (note zero-suppressed y axis) shows our results for temporal separations data in the two IC 59 data sets. Neither data set shows evidence for deviation from the expected flat distribution (illustrated using scrambled data sets), with ${\chi }^{2}$-derived p-values of $p=73 \%$ for IC 59-North and $p=53 \%$ for IC 59-South.

Zoom In Zoom Out Reset image size

Examining the angular and timing separations of the neutrinos and coincident photons at higher resolution thus provides no support for the presence of short-duration ($\delta t\lesssim 100$ s) ν + γ emitting cosmic sources as conceived in our original hypothesis. Since these are not seen, we move on to examine alternative models that might generate ν + γ spatio-temporal correlations in the data.

We carry out two tests for spatio-temporal correlations between the neutrino and gamma-ray data sets beyond our original ±100 s temporal acceptance window.

First, a correlation between neutrino and photon positions on the sky, without any temporal correlation (i.e., in steady state) could suppress λ values relative to the null hypothesis, due to the B_i gamma-ray background terms in our pseudo-likelihood (Equation (1)).

To test for positional correlation, we first construct a single Fermi background map covering the full energy range. We then measure the background value at the location of every IceCube neutrino in unscrambled data and compute the average photon background for the neutrino map. This average is then compared to the average backgrounds from each of the 10,000 scrambled data sets. The scrambled data sets give an average background of $(1.90\pm 0.015)\times {10}^{-2}$ photons deg⁻² m⁻² per 200 s for IC 59-North and $(2.40\pm 0.022)\times {10}^{-2}$ photons deg⁻² m⁻² per 200 s for IC 59-South. The observed average backgrounds (in the same units) from unscrambled data are $1.91\times {10}^{-2}$ (+0.58σ; $p\,=28.1 \%$) for IC 59-North and $2.44\times {10}^{-2}$ (+1.67σ; $p=4.7 \%$) for IC 59-South. These observed values are presented in the context of the distributions from scrambled data in Figure 9.

Zoom In Zoom Out Reset image size

This analysis is not an independent test for the presence of ν + γ sources, but rather, an attempt to identify an underlying reason for the trend in λ residuals seen in Section 4.1. Since the p-value for the separate analyses, as well as their combination ($p=7.0 \%$ by Fisher's method), are within a factor of two of the p-values from the corresponding λ distribution tests, this provides reason to interpret the latter result as (at least in part) due to the observed tendency of IC 59 neutrinos to land on systematically brighter regions of the gamma-ray sky. We reiterate that while this tendency is present in the data, it is not sufficiently strong (in a statistical sense) to support an evidence claim.

As an alternative approach, we check for correlated ν + γ variability on timescales beyond our predefined ±100 s temporal window but shorter than the full extent of the Fermi mission. To this end, for each neutrino in our scrambled IC 59 data sets, we use the total Fermi mission background map to calculate the number of photons expected to arrive within 5° of the neutrino position and ±50,000 s of the neutrino arrival time (excluding the ±100 s window used in the original analysis). We then count the number of photons arriving within this spatio-temporal window (again, summing results across our three Fermi energy bands). For each neutrino, the observed number of photons within the extended temporal window is expressed as a Poisson fluctuation on the number expected by normalizing against the full Fermi mission. We quantify the magnitude of this fluctuation as a p-value and find the equivalent number of σ for a Gaussian distribution, yielding a statistic that we call the local excursion E for that neutrino. The distribution of excursions from all neutrinos in unscrambled data can then be compared to expectations from scrambled data.

We perform the same two tests that we developed in our primary analysis above. First, we check for individual events that exhibit an unusually large excursion, exceeding either the 1 in 10 (${E}_{\times 10}$) or 1 in 100 (${E}_{\times 100}$) thresholds from scrambled data. Second, we compare the excursion distribution from unscrambled data to the null distribution from scrambled data using the Anderson–Darling k-sample test. Since this analysis involves two further independent tests of the two data sets, we add four trials for purposes of trial correction of any subsequent statistics, giving ${N}_{\mathrm{trials}}=12$.

Excursion thresholds for the two data sets are ${E}_{\times 10}=614$ and ${E}_{\times 100}=1285$ for IC 59-North, and ${E}_{\times 10}=333$ and ${E}_{\times 100}=1075$ for IC 59-South. As in our primary analysis, the highest-excursion events in the scrambled data are due to the two GRBs observed in the IC 59-North data. Excluding these GRBs would give excursion thresholds of ${E}_{\times 10}=575$ and ${E}_{\times 100}=1102$ for IC 59-North.

Analyzing the unscrambled IC 59 data sets reveals no excursions above the ${E}_{\times 10}$ threshold for either data set. Performing the Anderson–Darling test on the null and unscrambled distributions yields $p=55 \%$ for IC 59-North and $p=62 \%$ for IC 59-South. We therefore see no evidence for spatio-temporal correlation of the neutrinos and Fermi gamma rays on the ∼0.5 day timescale that we probed.

We conclude that the observed tendency of IC 59 neutrinos to arrive from brighter portions of the Fermi gamma-ray sky, while potentially due to a statistical fluctuation (single-trial $p=7.0 \%$ for the two hemispheres combined), is both intriguing in its own right and likely explains the systematic trends in λ residuals against scrambled data sets observed for both hemispheres in our original analysis (single-trial $p=3.9 \%$ for the two hemispheres combined). While this p-value cannot support an evidence or discovery claim in the context of our multistage analysis, it nonetheless points the way to interesting future analyses that could make use of eight further years of IceCube data from the 79-string and full-strength (86-string) arrays.

In particular, we note that it is a low-level (single neutrino) correlation between the neutrino and gamma-ray skies that has prompted current interest in the blazar TXS 0506 + 056 and its possible neutrino IceCube-170922A (Tanaka et al. 2017).