Searches for IceCube neutrinos from blazar flares using correlations with gamma-ray lightcurves

Blazars have long been considered as accelerator candidates for cosmic rays. In such a scenario, hadronic interactions in the jet would produce neutrinos and gamma rays. Correlating the astrophysical neutrinos detected by IceCube with the gamma-ray emission from blazars could therefore help elucidate the origin of cosmic rays. In our method we focus on periods where blazars show an enhanced gamma-ray flux, as measured by Fermi-LAT, thereby reducing the background of the search. We present results for TXS 0506+056, using nearly 10 years of IceCube data and discuss them in the context of other recent analyses on this source. In addition, we give an outlook on applying this method in a stacked search for the combined emission from a selection of variable Fermi blazars.


Introduction
Blazars are a promising class of potential cosmic-ray accelerators and neutrino sources. They have been observed, e.g. by Fermi-LAT [1], as highly variable sources of high-energy gamma rays, which may be produced together with neutrinos through the p-p and p-γ interactions of cosmic rays [2]. Neutrinos would easily escape and provide a "smoking gun" for the presence of these hadronic mechanisms. IceCube is a neutrino telescope deployed in a cubic kilometre of glacial ice at the geographic South Pole [3]. In the case of a charged-current interaction of a muon neutrino its direction can be reconstructed with a resolution better than 1 • for E ν > 1 TeV [4]. A diffuse flux of astrophysical ν µ has been measured [5], but no individual point source discovered [4]. The coincidence between a single Extremely High Energy (EHE) neutrino registered by IceCube on September 22nd 2019 and a blazar flare from TXS 0506+056 [6] could be the first evidence for a astrophysical neutrino source. As can be seen in figure 1, the 2017 flare began several months prior to the neutrino event and the same blazar was also active in the years before. This motivates a search for additional neutrinos which correlate with the lightcurve, as discussed hereafter. Similar analyses have been performed on other blazars individually [7], however the method had not yet been used in the kind of stacking analyses which are the most powerful probes of the contribution of blazars to the astrophysical neutrino flux [8]. In this contribution we present such a combined analysis which could provide stronger limits, or a discovery.

Method
Lightcurve correlation We use the maximum likelihood method described in [7]. From the lightcurve LC(t) it derives a signal time probability distribution function (PDF), which gets combined with the point spread function (PSF) and energy PDF p(E) to define the signal term S in the likelihood: x, E, and t are respectively the spatial coordinates, energy, and arrival time of a given neutrino event. The fitted parameters in the analysis are the number of signal events (n S ), the spectral index γ of the neutrino energy spectrum, and the photon flux threshold Φ 0 . For Φ 0 = 0, this is just the normalized lightcurve i.e. all gamma rays and neutrinos are produced together. For fitted values of Φ 0 > 0, only the excess above Φ 0 contributes to neutrino signals. Flare stacking We give sources thresholds by combining a common threshold parameter τ with lightcurve-specific metrics q, for the steady quiescent flux, and σ q , for its r.m.s., as defined in [9]: To stack, the signal terms of multiple sources are summed as in S (x, E, t; γ, τ) = k∈sources w (k) (γ, τ)S k (x, E, t; γ, Φ (k) 0 (τ)). The relative source weight w (k) includes one factor from the exposure of the detector to a spectrum of the given power-law index γ at the source's declination [8]. We also include factors from three alternative weighting schemes, which depend on τ: 1. Where neutrinos and gamma rays are produced by pion decays, one obtains L ν ∝ L γ [2]. Assuming that the relative effects of gamma-ray reprocessing and the finite energy bands covered by Fermi-LAT and IceCube respectively are constant, this leads to a weighting proportional to the energy flux measured in Fermi-LAT. 2. From p − γ interactions in a blazar where the amount of injected protons are themselves proportional to the photon luminosity, one obtains L ν ∝ L 2 γ as summarised in [10]. In this case, the luminosity distance d L is needed to calculate the luminosity explicitly.
3. An equal weighting, which does not reflect a physical model itself, but is a generic case that can provide sensitivity to a true physical scenario very misaligned with the other weighting schemes.
Calculating limits and sensitivities In order to estimate the sensitivity of the analysis as well as the significance (p-value) of the result, simulated data sets are needed. We follow a similar approach as [4]. Real data events are used as background events, by re-assigning their arrival times to random times selected from periods when the detector was operating, then re-computing equatorial coordinates. Signal events are added from detector simulation data sets, with energies according to the assumed energy spectrum and arrival times from the signal time distribution. Also for that, only the detector operation periods are considered.

Gamma-ray and neutrino data samples
Fermi-LAT lightcurves We use lightcurves based on Fermi-LAT PASS8 data 1 . Two strategies are used to avoid contamination by nearby sources. First, the photon energy threshold is raised, improving the angular resolution. Then, a likelihood fit on each time bin optimizes not only the normalisation of the observed source, but also of its neighbours. For the blazar flare stacking we use lightcurves used in [6] with an energy threshold of 1 GeV and a 28-day binning. For the study of TXS 0506+056, the energy threshold is 700 MeV and the binning 7 days. In both cases we smooth the statistical fluctuations with Bayesian Blocks [11]. Source selection For the blazar flare stacking we select from 2254 extragalactic Fermi sources which were used in [6] where the lightcurve contains a minimum number of goodquality fits and has at least a minimum amount of variability. This leaves 114 Flat Spectrum Radio Quasars (FSRQs) and 65 BL Lacs, as classified by the Fermi-LAT 3FGL catalog [12]. Since these two source classes have different luminosity ranges, we separate them in the analysis. When using the L 2 γ weighting scheme (see sec. 2), we remove 19 sources for which no redshift could be found in TeVCat 2 , 3FGL, or NED 3 . Neutrino data We use the same data as the searches for neutrinos from TXS 0506+056 [13].

TXS 0506+056 lightcurve correlation analysis
To separate the outcome of this analysis from the observation of IceCube-170922A that triggered it, this event is removed from the data sample. Maximizing the likelihood described in sec. 2 results in the left column of table 1. Using the method of sec. 2 the resulting p-value is 12.6%. As a separate check, the IceCube-170922A event is re-inserted into the data and the fit performed again, with results as in the right column of table 1 and a p-value of 0.065% 4 . Figure 2 shows limits on the number of neutrino events with a spectral index γ = 2 and a temporal distribution described by applying a range of thresholds. In the most conservative case of Φ 0 = 0, we find a constraint of <11 events (apart from IceCube-170922A). For higher thresholds, signal events must concentrate during the shorter period of higher flux and therefore limits improve.
Blazar flare stacking In this analysis TXS 0506+056 is removed since this source was similarly analysed on its own (see sec. 4). Table 2 shows computed sensitivities at 90% confidence level in terms of number of signal events. The 64 brightest FSRQs are taken as a proxy for the complete list due to limited computing resources, but all 114 will be used for the eventual results. A benchmark γ = 2 and conservative τ = −∞ are used.  Table 2. Sensitivities for the blazar flare stacking analysis at 90% C.L. in terms of number of signal events from the hypothesis γ = 2, τ = −∞.

Discussion
This analysis looking for neutrino events correlated to the gamma-ray lightcurve of TXS 0506+056 is sensitive to the EHE neutrino event IceCube-170922A which prompted its study. However, when removing this event, no additional significant signal is found. It is important to note that these limits do not apply to neutrinos that do not correlate with the gamma-ray lightcurve, such as the excess in 2014 [13], indicated in figure 1. It remains of interest how these three observations will be unified. The stacking analysis described here provides a powerful means of identifying the contribution of neutrinos from blazar flares to the total cosmic neutrino flux. This may help resolve the role of blazars as cosmic-ray sources.