SEARCHES FOR TIME-DEPENDENT NEUTRINO SOURCES WITH ICECUBE DATA FROM 2008 TO 2012

The cosmic-ray spectrum spans 10 decades in energy up to 10¹¹ GeV per particle, and despite being extensively studied for many years, the origin and the acceleration mechanism remain uncertain. Cosmic rays consist of hadrons, mainly protons, and in part also ionized nuclei with an energy-dependent composition. As a result of interactions of cosmic rays with matter and ambient photons close to the acceleration sites, pions are produced, and decays of the charged pions and their daughter muons then produce neutrinos. Such astrophysical neutrinos are a unique and valuable messenger in astroparticle physics because of their properties and the specifics of their production mechanisms. In particular, they carry information about the origin and spectrum of cosmic rays. In contrast to the cosmic rays, neutrinos are not deflected in magnetic fields, nor do they interact on the way to Earth. As a consequence, their trajectories point back to their origin. A detection of neutrinos from a given site would therefore stand as proof of accelerated hadrons, identifying it as a source of cosmic rays.

This paper updates and expands the IceCube time-dependent searches and results for flaring sources (Abbasi et al. 2012a, 2012b; Aartsen et al. 2013a) and periodic sources (Abbasi et al. 2012c). Time-dependent astrophysical neutrino signals can be better observed by using event times in addition to the direction and energy used in standard IceCube point-source searches, because the additional information improves the rejection of the atmospheric muon and neutrino events that form the dominant background.

If several events are coming from the same astrophysical point source, they should be spatially concentrated around the emitting source, and they should have a harder spectrum than muons and neutrinos produced in atmospheric showers. Assuming a Fermi acceleration model (Bell 1978; Schlickeiser 1989), a differential neutrino spectrum close to ${E}^{-2}$ will be produced. Additional effects like the acceleration of muons in the cosmic-ray sources (Klein et al. 2013; Winter et al. 2014) can modify the spectrum. In case of sufficiently high acceleration gradients (above 160 keV m⁻¹), often required for extremly short flares, the daughter muons from charged pion decays can be accelerated before they decay. Thus, the energy of the neutrinos from the muon decays will be enhanced. The impact on the overall neutrino spectra strongly depends on various properties of the source, e.g., the magnetic field strength and the amount of matter in the acceleration region. Conventional atmospheric neutrinos at energies above 100 GeV have a much softer spectrum that asymptotically approaches a spectrum one power steeper than the primary spectrum, first in the vertical direction and at higher energies also for larger zenith angles (Gaisser 2005; Aartsen et al. 2014).

Time-dependent searches can be performed in an "untriggered" way, meaning that the correlation of event times is investigated and no additional information from independent observations is used. The full parameter space of time, energy, and direction of measured events is scanned looking for clusters in time and space of high-energy events among the background of atmospheric events, hereby called the "All-Sky Time Scan." This search is the most generic one, but it is subject to a large trial factor that penalizes the significance of the signal with respect to the background. Hence, in addition, more specific searches are carried out using a multi-messenger approach. The basic assumption is that neutrinos and γ-rays are correlated because they have a common origin in the astrophysical sources that accelerate the cosmic rays. The "Search for Triggered Multi-Messenger Flares" makes use of light curves measured by γ-ray experiments. This analysis is triggered by multiwavelength measurements (alerts for flares from TeV and X-ray experiments), as well as the Fermi-LAT light curves, which provide continuous monitoring of selected sources or of flares above a certain photon flux level.⁵¹

The paper outline is as follows. In Section 2 the possible sources of astrophysical neutrinos are briefly discussed, and in Section 3 the IceCube detector and the details of the data samples used are introduced. Section 4 gives a detailed description of the general likelihood method applied while the specifics, which are different for each search, are given at the beginning of the corresponding sections.

In the following sections, the results for the various searches are presented. First, in Section 5 the untriggered "All-Sky Time Scan" looking for any cluster of high-energy neutrino events from any direction in the sky is presented. For this search the data taken by the 59-string configuration of IceCube (IC-59), the 79-string configuration (IC-79), and the first year of the full IceCube data taking with 86 strings (IC-86I) were used. For each of the three years time-dependent skymap scans were performed.

Then in Section 6 the results of the "Search for Triggered Multi-Messenger Flares" are discussed. The data taken with the IC-59, IC-79, and IC-86I configurations of IceCube were combined. Hence, this represents a long-term study for flaring objects. In Section 7 the results of the "Search for Triggered Flares with Sporadic Coverage" are presented. These were triggered by higher-energy γ-ray observations in the TeV range.

In Section 8 the "Search for Periodic Neutrino Emissions from Binary Systems" is presented. The search was performed in the phase domain. For this search data from four IceCube data-taking seasons, from the IC-40 to the IC-86I configuration, were used.

Finally, in Section 9 we conclude with a brief outlook.

The focus of this paper is on active galactic nuclei (AGNs), and particularly blazars, i.e., AGNs with jets pointing toward us, which are interesting due to their variability and the high power emitted during bursts (Begelman et al. 1984; Antonucci 1993; Urry & Padovani 1995). The analyses described in the following sections are sensitive also to short flares (of the duration of a few seconds) and hence, to some extent, to γ-ray bursts. More sensitive searches for neutrinos from γ-ray bursts are presented elsewhere (Abbasi et al. 2011a, 2012d). Blazars exhibit sudden sequences of multiple flares that may last from minutes to months and are observed in various wavelengths from radio to γ-rays. Correlations between various bands have been observed in numerous multiwavelength campaigns, particularly between X-rays and γ TeV emissions. A well-studied case is that of the close-by TeV blazar Mrk 421 (see, e.g., Acciari et al. 2011; Aleksic et al. 2012b). In other instances optical flares have triggered TeV flare observations (Reinthal et al. 2012).

The spectral energy distribution of AGNs is characterized by two broad peaks. The lower-energy one is believed to originate from synchrotron emission of the charged particles in the jet. In leptonic models the higher-energy one is generally explained by inverse Compton scattering of either the synchrotron seed photons (synchrotron self-Compton [SSC]; see, e.g., Maraschi et al. 1992) or external seed photons (external Compton [EC]; Dermer & Schlickeiser 1993; Ghisellini et al. 2005) by the electrons and positrons in the jet.

In the simplest case, both SSC and EC mechanisms predict that flaring at TeV energy should be accompanied by a simultaneous flaring in the synchrotron peak, and so a connection between bands from optical to γ-ray is expected. If the synchrotron peak is located far from the optical, as in the case of high-frequency peaked blazars, then the synchrotron flares should be visible at other wavelengths, usually X-rays (Acciari et al. 2011; Aleksic et al. 2012b).

Alternatively, hadronic models suggest the production of neutrinos and γ-rays from pion decays (Mannheim 1993; Atoyan et al. 2002; Mücke et al. 2003). In these models, the high-energy peak is due to proton synchrotron emission or decay of neutral pions formed in cascades by the interaction of the high-energy proton beam with the radiation or gas clouds surrounding the source (Atoyan et al. 2002). In this scenario, a strong correlation between the γ-ray and the neutrino fluxes is expected. Observations of neutrinos would clearly distinguish between leptonic and hadronic models.

For some observations it has been claimed that hadronic processes could explain flares better than leptonic processes (Böttcher et al. 2009). Orphan flares, i.e., a TeV flare without a lower-energy counterpart, challenge leptonic models (Halzen & Hooper 2005; Reimer et al. 2005; Sahu et al. 2013). Non-observation of significant X-ray activity could naturally be interpreted as due to the suppression of electron acceleration and inverse Compton scattering and dominance of very high energy (VHE) γ-ray production from meson decays in hadronic models.

In addition to flares, another possible time-dependent signature is a periodicity in the neutrino emission. In the case of binary systems, a modulation of the neutrino emission could occur due to the relative geometry of the emitted beam with respect to the large-mass star and the observer. A particularly interesting case is that of microquasars, which are radio-jet X-ray binaries that can include either a neutron star or a black hole. A similar modulation is observed for three binary systems in TeV γ-rays: LS 5039 (Aharonian et al. 2006a), LS I+61 303 (Smith et al. 2013), and HESS J0632+057 (Aliu et al. 2013). These three binary systems are found to be periodic at GeV and TeV energies, although the emission in the two bands seems anticorrelated (Hadasch et al. 2012). This anticorrelation is in fact a generic feature of models in which the GeV (TeV) emission is enhanced (reduced) when the highly relativistic electrons are moving in the direction of the observer and encounter the seed photons head-on. Neutrino emission from microquasars has been described in various papers (Distefano et al. 2002; Romero et al. 2005; Aharonian et al. 2006b; Christiansen et al. 2006; Torres & Halzen 2007). It can be assumed that a neutrino signal would be correlated to the periodic γ-ray emission from the binary system due to the system's rotation, but it is not clear in which phase with respect to the γ-rays the neutrinos are to be expected. Hence, the values for the period are taken from multiwavelength observations, and the phase is fitted as a free parameter.

The IceCube observatory is a Cerenkov detector searching for high-energy neutrinos. It is an array of 5160 digital optical modules (DOMs) deployed deep in the Antarctic ice near the South Pole. The purpose of the detector is to allow observations of neutrinos of astrophysical origin and atmospheric muons and neutrinos induced by cosmic rays at energies around and above the knee ($\sim 3\times \ {10}^{15}\;\mathrm{eV}$). Cerenkov light produced by secondary leptons from neutrino interactions in the vicinity of the detector is used to indirectly detect these neutrinos. Only charged-current muon neutrino events were considered for the studies presented here, because of the large range of the secondary muons, which provides good angular resolution and a high rate by including events that start outside the detector. The pointing information relies on the secondary muon direction, which at energies above TeV differs from the original neutrino direction by less than the angular resolution of the detector (Aartsen et al. 2013b).

The DOMs are spherical, pressure-resistant glass housings, each containing a Hamamatsu photomultiplier tube of 25 cm diameter and the associated electronics necessary for waveform digitization (Abbasi et al. 2009, 2010). These DOMs are connected together in 86 vertical strings of 60 DOMs each, and the strings are lowered into 60 cm wide holes drilled into the ice using hot water to instrument a kilometer-cubed volume in a depth range from 1.5 down to 2.5 km. At the center of the detector the strings are placed closer to each other, forming a sub-array called DeepCore. This sub-array enhances the sensitivity of the detector to neutrinos of energies below the standard IceCube threshold of 100 GeV (Abbasi et al. 2012e).

The IceCube detector has been in operation even before full completion in 2010 December. The partial detector layouts (resp. data samples) used for this paper are labeled IC-40, IC-59, IC-79, and IC-86I, where the number corresponds to the number of operational stings. The data samples used in the studies presented here are described in detail in previous papers reporting time-integrated searches (Abbasi et al. 2011b; Aartsen et al. 2013b, 2013c).

The analyses in the present paper used identical detector simulation and corresponding estimates of neutrino effective area and point-spread function based on Monte Carlo studies as for the time-integrated analyses (Aartsen et al. 2013b). Table 1 summarizes the four data samples consisting of both up- and down-going muons, which have different origins; while the up-going muons are mostly from interactions of atmospheric neutrinos that have passed through the Earth, the down-going muons are mainly from meson decays in atmospheric showers caused by cosmic rays.

The unbinned time-dependent likelihood method that was used for the analyses in this paper has been applied in previous analyses (Braun et al. 2008) and (Braun et al. 2010).

The likelihood function is defined as

$$\begin{eqnarray}{\mathcal{L}}(\gamma ,{n}_{{\rm{s}}},\ldots ) & = & \underset{j}{\displaystyle \prod }{{\mathcal{L}}}^{j}(\gamma ,{n}_{{\rm{s}}}^{j},\ldots )\\ & = & \underset{j}{\displaystyle \prod }\;\underset{i\in j}{\displaystyle \prod }\left[\displaystyle \frac{{n}_{{\rm{s}}}^{j}}{{N}^{j}}{{\mathcal{S}}}_{i}^{j}+\left(1-\displaystyle \frac{{n}_{{\rm{s}}}^{j}}{{N}^{j}}\right){{\mathcal{B}}}_{i}^{j}\right],\end{eqnarray} \tag{ 1 }$$

where $i\in j$ indicates that the ith event belongs to the jth sample (IC-40, IC-59, IC-79, or IC-86I), N^j is the total number of events in the jth sample, ${{\mathcal{S}}}_{i}^{j}$ and ${{\mathcal{B}}}_{i}^{j}$ are the signal and background probability density function (pdf), respectively, and γ is the spectral index of the differential spectrum ${dN}/{dE}\propto {E}^{-\gamma }$ and is assumed to be the same for all the data samples. The total number of signal events ${n}_{{\rm{s}}}$ is the sum of the signal events coming from each sample ${n}_{{\rm{s}}}^{j}$. The fractional contribution of each ${n}_{{\rm{s}}}^{j}$ to the total ${n}_{{\rm{s}}}$ is fixed by the relative neutrino effective areas of each of the configurations (determined by detector simulation) and varies depending on the spectral index and the declination. The likelihood function in Equation (1) is a function of the global common parameters γ and ${n}_{{\rm{s}}}$ and the parameters specific for each time-dependent search.

The signal pdf ${{\mathcal{S}}}_{i}^{j}$ is given by

$$\begin{eqnarray}&&{{\mathcal{S}}}_{i}^{j}={P}_{i}^{\mathrm{signal},j}(| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{{\rm{s}}}| ,{\sigma }_{i})\cdot {{\mathcal{E}}}_{i}^{\mathrm{signal},j}({E}_{i},{\delta }_{i},\gamma )\cdot {{\mathcal{T}}}_{i}^{\mathrm{signal},j},\end{eqnarray} \tag{ 2 }$$

where ${{\boldsymbol{x}}}_{{\rm{s}}}$ is the source direction, ${{\boldsymbol{x}}}_{i}$ is the reconstructed direction of the event, ${\sigma }_{i}$ is the angular error estimate for the reconstruction, and E_i and ${\delta }_{i}$ are the reconstructed energy and declination, respectively. The term ${P}_{i}^{\mathrm{signal},j}$ is the spatial pdf (the point-spread function), and ${{\mathcal{E}}}_{i}^{\mathrm{signal},j}$ is the energy pdf. These two terms are the same as for the IceCube time-integrated searches and can be found in Aartsen et al. (2013c). The time-dependent signal pdf ${{\mathcal{T}}}_{i}^{\mathrm{signal},j}$ is specific for each of the different signal hypotheses and will be described in the corresponding sections.

Since the data are background dominated, the background pdf ${{\mathcal{B}}}_{i}^{j}$ can be estimated from the data themselves, and it has the form

$$\begin{eqnarray}&&{{\mathcal{B}}}_{i}^{j}={P}_{i}^{\mathrm{bkg},j}({\delta }_{i})\cdot {{\mathcal{E}}}_{i}^{\mathrm{bkg},j}({E}_{i},{\delta }_{i})\cdot {{\mathcal{T}}}_{i}^{\mathrm{bkg},j}.\end{eqnarray} \tag{ 3 }$$

For sufficiently long timescales, the spatial pdf for the background ${P}_{i}^{\mathrm{bkg},j}({\delta }_{i})$ is uniform in right ascension as the Earth rotation averages over detector effects and follows the distribution of the data with declination. For short timescales (less than a day) the spatial pdf is no longer uniform in right ascension and shows a pattern when the neutrino directions align with the directions of geometrical symmetries of the non-uniform detector array. The effect disappears as the timescale grows large enough to allow the Earth rotation to average out the differences. The energy density function ${{\mathcal{E}}}_{i}^{\mathrm{bkg},j}$ can be found in Aartsen et al. (2013c). The time probability ${{\mathcal{T}}}_{i}^{\mathrm{bkg},j}$ is taken to be flat since it was demonstrated in Abbasi et al. (2012a) that the seasonal modulations of background atmospheric muons and neutrinos are negligible compared to possible measurable signals.

For all the searches presented here the test statistic TS is defined as the maximum likelihood ratio:

$$\begin{eqnarray}&&{TS}=-2\mathrm{log}\left[\displaystyle \frac{{\mathcal{L}}({n}_{{\rm{s}}}=0)}{{\mathcal{L}}({\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},\ldots )}\right],\end{eqnarray} \tag{ 4 }$$

where ${\mathcal{L}}({n}_{{\rm{s}}}=0)$ is the likelihood of the background-only (null) hypothesis and ${\mathcal{L}}({\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},\ldots )$ is the likelihood evaluated with the best-fit values for the signal-plus-background hypothesis. The parameters denoted with carets are the best-fit values of ${n}_{{\rm{s}}}$ and ${\gamma }_{{\rm{s}}}$, which are present in all searches,⁵² while the parameters specific to each time-dependent search (represented by the ellipsis above) will be discussed in each search section. Estimates for the values of these parameters can be derived from the maximization of the likelihood function defined in Equation (1).

The test statistic TS as defined in Equation (4) for the background-only data is expected to follow a χ² distribution with number of degrees of freedom reflecting the number of fitted parameters. This behavior of the TS can be used to estimate the pre-trial p value as the probability of the TS value from the ${\chi }^{2}$ distribution.

Because our samples are background dominated, in order to translate the pre-trial p values into post-trial p values, the time-scrambled IceCube data from the same period can be used to generate background samples. The time-scrambling procedure assigns to the events a random time within the period while keeping all other event properties (energy, local reconstructed coordinates, etc.) unchanged. In this way many different samples describing the background are obtained. About 1000 of these sets were generated, and for each of them a p value was obtained. Then the distributions of these p values for the time-scrambled data were compared to the p value obtained for real data. From this the probability is derived that a random fluctuation will result in deviation at least as significant as our result for real data, and this probability is called the post-trial p value. The obtained post-trial p value is stable with respect to deviations of background-only TS from the assumed ${\chi }^{2}$ distribution and accounts for trial factors related to looking at many positions at the sky.

In analyzing the results, it is convenient to characterize events by their time-integrated weight w_i, defined as the ratio of the signal and background pdfs from Equations (2) and (3) without the time pdf term:

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{{P}_{i}^{\mathrm{signal},j}(| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{{\rm{s}}}| ,{\sigma }_{i})\cdot {{\mathcal{E}}}_{i}^{\mathrm{signal},j}({E}_{i},{\delta }_{i},\gamma )}{{P}_{i}^{\mathrm{bkg},j}({\delta }_{i})\cdot {{\mathcal{E}}}_{i}^{\mathrm{bkg},j}({E}_{i},{\delta }_{i})},\end{eqnarray} \tag{ 5 }$$

which will be later used to visualize the result in a way helpful for understanding whether the significance comes from the spatial or the time clustering.

For the time-dependent searches⁵³ it is natural to express the discovery potentials, sensitivities, and upper limits in terms of fluence, defined for an ${E}^{-2}$ spectrum as

$$\begin{eqnarray}f & = & {\displaystyle \int }_{{t}_{\mathrm{min}}}^{{t}_{\mathrm{max}}}{dt}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}\times E\displaystyle \frac{{dN}}{{dE}}\\ & = & {\rm{\Delta }}t{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}\times E\displaystyle \frac{{{\rm{\Phi }}}_{0}}{{E}^{2}}\\ & = & {\rm{\Delta }}t{{\rm{\Phi }}}_{0}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{1}{E}{dE}\end{eqnarray} \tag{ 6 }$$

where ${{\rm{\Phi }}}_{0}$ is the normalization of an average flux during the emission. The emission duration ${\rm{\Delta }}t$ is defined by the time pdf. The integration limits ${E}_{\mathrm{min}}$ and ${E}_{\mathrm{max}}$ are set to the 5% and 95% energy percentile, respectively, of the event sample for the given declination band.

This is the most generic time-dependent search among the ones presented in this paper. It is optimized to look for neutrino emission from a point-like source with limited duration. Because it aims for one-time flares, it does not benefit from adding multiple data sets together. Here the results for the IC-59, IC-79, and IC-86I IceCube data are presented separately. The results of a similar search for the IC-40 configuration can be found in Abbasi et al. (2012a).

5.1. Method

The "All-Sky Time Scan" is an untriggered search since it is performed only using the neutrino data themselves. For this search the whole sky⁵⁴ was divided into a grid of $0\buildrel{\circ}\over{.} 1\times 0\buildrel{\circ}\over{.} 1$, much finer than the angular resolution of the detector, and the likelihood was maximized at each grid point. The expression for the likelihood is given by Equation (1) (with each detector configuration analyzed separately). In Equation (2) the signal pdf ${{\mathcal{T}}}_{i}^{\mathrm{signal}}$ was chosen of the following form:

$$\begin{eqnarray}&&{{\mathcal{T}}}_{i}^{\mathrm{signal}}=\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{T}}\mathrm{exp}\left(-\displaystyle \frac{{({t}_{i}-{T}_{0})}^{2}}{2{\sigma }_{T}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

where t_i is the arrival time of the ith event, and T₀ and ${\sigma }_{T}$ are the mean and the width of a Gaussian time structure, respectively. The Gaussian function was used as a smooth, general parametrization of a limited duration increase in the emission of a source. The best-fit values of these parameters were obtained maximizing the likelihood function described in Section 4.

Since within the limited livetime of a given sample it is possible to accommodate a larger number of flares as their duration decreases, an undesired bias toward finding a short flare is introduced. This could also be interpreted as a hidden extra trial factor affecting short flares. To account for it, the test statistic TS from Equation (4) was modified. As explained in Braun et al. (2010), an additional marginalization term $T/\sqrt{2\pi }\sigma$ that penalizes short flares compared to long ones was introduced. The test statistic including this modification is thus

$$\begin{eqnarray}&&{TS}=-2\mathrm{log}\left[\displaystyle \frac{T}{\sqrt{2\pi }{\hat{\sigma }}_{T}}\times \displaystyle \frac{{\mathcal{L}}({n}_{{\rm{s}}}=0)}{{\mathcal{L}}({\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{\sigma }}_{T},{\hat{T}}_{0})}\right],\end{eqnarray} \tag{ 8 }$$

where ${\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{\sigma }}_{T},{\hat{T}}_{0}$ are the best-fit values and T is the total livetime of the data-taking period (either IC-59, IC-79, or IC-86I). For details on the numerical maximization procedure see Abbasi et al. (2012a).

The expected performance of this approach in terms of discovery potential and sensitivity is shown in Figure 1. The figure compares this time-dependent search with the standard time-integrated all-sky search for given declination. The "All-Sky Time Scan" search is better for flares with ${\sigma }_{T}$ shorter than 100 days in terms of discovery potential and for flares with ${\sigma }_{T}$ shorter than 6 hr in terms of sensitivity. As the flare duration ${\sigma }_{T}$ gets shorter, the sensitivity levels out at around three events. For the calculation of the discovery potential at least two events are required in order to identify a flare. To generate a sample according to a Poissonian distribution with 50% of the cases having two or more events, the Poissonian mean has to be equal to 1.68. Therefore, for the shortest timescales the discovery potential will asymptotically approach this value, causing it to drop below the sensitivity.

Zoom In Zoom Out Reset image size

5.2. Results

The maximization of the likelihood at each grid point in the sky results in a map of TS values or, equivalently, a map of pre-trial p values that serve as an estimate of the local significance of the best-fit parameters. To address the question of whether any excess in the sky is significant, a correction for the trial factor involved in searching the entire sky had to be used. This was done as described at the end of Section 4 by repeating the analysis on time-scrambled data.

5.2.1. IC-59 Results

Figure 2 shows the IC-59 skymap of pre-trial p values for the all-sky search. The most significant point in the IC-59 data was found at (R.A., decl.) = (2135, −025). The peak occurred on MJD 55,259 (2010 March 4) and had a width parameter ${\hat{\sigma }}_{T}$ of 5.5 days, a soft spectral index of $\hat{\gamma }=3.9$, and 14.5 fitted signal events. The pre-trial p value was $2.04\times {10}^{-7}$; a value at least as significant as this was found somewhere in the sky in 14 out of 1000 scrambled maps. Thus, the post-trial p value was 1.4%, which was low but not significant evidence of an actual flare. When this analysis was repeated on the data for the following years (these results are presented in the following sections), the IC-59 hot spot significance decreased, and it was not seen with high significance any more. Figure 3 shows the time-integrated event weights w_i at the position of maximum significance plotted throughout the year; a clustering near the time of the best-fit ${\hat{T}}_{0}$ is clearly visible. When a bin of radius $2^\circ$ and 13 days in time (the FWHM of the Gaussian) centered on the peak is considered, 13 events are found compared to an expected background of 1.7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2.2. IC-79 Results

Figure 4 shows the IC-79 skymap of pre-trial p values. The most significant deviation from the background-only hypothesis was found at (R.A., decl.) = (34345, −3165). The mean of the fitted Gaussian flare was at MJD 55,466 (2010 September 27) with a width parameter ${\hat{\sigma }}_{T}$ of 1.8 days, a soft spectral index of $\hat{\gamma }=3.95$, and 7.2 fitted signal events. The large blue spot in the upper right quadrant in Figure 4 was caused by two events arriving very close in time, the consequence of which was an increased significance over a wider area to which those two events contribute.

Zoom In Zoom Out Reset image size

The pre-trial p value obtained for the IC-79 hot spot was $1.07\times {10}^{-5}$. The post-trial p value was 66%, consistent with a background-only hypothesis.

Figure 5 shows the time-integrated event weights at the position of maximum significance plotted throughout the year with a hint of clustering recognizable at the fitted time.

Zoom In Zoom Out Reset image size

5.2.3. IC-86I Results

Figure 6 shows the skymap of the pre-trail p values obtained for the IC-86I data; the most significant point was at (R.A., decl.) = (23595, 4295), and the fitted Gaussian parameters are MJD 55,882 for the mean and 7.57 days for the width parameter ${\hat{\sigma }}_{T}$. A spectral index of $\hat{\gamma }=2.0$ and 13.1 signal events were fitted. A pre-trial p value of $1.06\times {10}^{-5}$ translates into a post-trial p value of 63%. The time-integrated event weights for the region around the most significant point through the year are shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This search targets a set of astronomical objects that were observed to be in a flaring state by Fermi-LAT⁵⁵ during the time period analyzed in this study. The tested hypothesis is that the neutrino emission follows the intensity of the photon emission. If the photon emission has a hadronic origin, then the accelerated hadrons (protons or nuclei) interact with matter and produce neutral and charged pions, which produce γ-rays and neutrinos, respectively, when they decay.

Since Fermi-LAT provides continuous monitoring and the data are publicly available, it is possible to make use of the measured light curves. A continuation of the search made during the previous period of IceCube data taking (Abbasi et al. 2012a; Aartsen et al. 2013a) is presented here, and it was extended including the whole IC-59, IC-79, and IC-86I combined data samples. Also a more advanced denoising method was implemented, which is described below, to better reconstruct the Fermi-LAT light curves.

6.1. Method

This search was performed only for selected objects in the sky. The criteria for selecting these objects (FSRQs, BL Lacs, etc.) were based on the Fermi-LAT photometric measurements. The Fermi-LAT monitored source list⁵⁶ was taken as a starting point. The first step was to retrieve the raw light curves from the Fermi Public Release data, using the analysis tools made available by the Fermi-LAT Collaboration. For each source the Fermi Science Tools v9r31p1 package⁵⁷ was used to select photons within 2° of the source and to calculate the exposure. Photon events with zenith angles greater than 105° were excluded to avoid contamination due to the Earth's albedo. For each source light curves with 1 day binning were obtained.

The denoised light curve was used as time-dependent signal pdf, ${{\mathcal{T}}}_{i}^{\mathrm{signal},j}$. The "Bayesian blocks" method (Resconi et al. 2009; Scargle et al. 2012) was applied to denoise the light curves, implemented for the purpose of this analysis in the version described in Scargle et al. (2012). A simplified explanation of the method is that it splits the time axis of the light curve into blocks for which the flux variations are assumed to be compatible with Poisson distributions with a constant mean within each block. The criterion for deciding whether to split a section into two blocks or not is based on comparing how well the variations follow a single Poisson distribution or two Poisson distributions with different means.

In order to optimize the performance of the Bayesian blocks method on IceCube data, the optimal value of the method parameter ${F}_{{\rm{B}}}$ had to be found. This parameter affects the method behavior in the following way. If the log-likelihood of two Poisson distributions is larger than the log-likelihood value for a single Poisson distribution by at least ${F}_{{\rm{B}}}$, the split is made. Too small values of ${F}_{{\rm{B}}}$ will cause denoised light curves to follow almost every point in the light curve, while for too large values of ${F}_{{\rm{B}}}$ the denoised result will ignore important structures of the light curve. In Figure 8 two examples of the denoised light curves are shown for ${F}_{{\rm{B}}}$ being outside of the optimal range. To determine the value of the parameter ${F}_{{\rm{B}}}$ to be used in this analysis, a series of tests were performed to evaluate the performance of the Bayesian blocks method as a function of ${F}_{{\rm{B}}}$.

Zoom In Zoom Out Reset image size

Realistic Fermi-LAT light curves were simulated in order to optimize the ${F}_{{\rm{B}}}$ parameter, using Fermi-LAT exposure data and assuming constant flux levels with Poissonian fluctuations. Folding in the Poisson-distributed flux with the exposure gave a model of the background fluctuations including the correct statistical errors. On top of this background Gaussian-shaped flares with random mean in time were injected. The background level was set to $0.5\times {10}^{-6}\;\mathrm{photons}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ and was not varied. Instead, the properties of the Gaussian flares were varied, i.e., strength of the flare, the time, and the duration. The amplitude values were taken to be either equal to the background level or twice the background level. The tested widths were 2 days, 5 days, and 10 days. In addition, combinations of two injected flares were also tested.

For each of the tested combinations of a flare with some amplitude and duration, hundreds of random instances of light curves with injected flares and background were produced. For each of these the light curve was denoised, scanning over different values of ${F}_{{\rm{B}}}$, and two quantities were calculated as a function of ${F}_{{\rm{B}}}$: the rate of finding a fake flare, i.e., a flare at a position where none was injected; and the rate of finding the injected flare. For this purpose "successful flare finding" was defined as the denoised light curve reaching above a three-standard-deviation upward fluctuation of the background. In Figure 9 an example is shown of two flares of different duration being injected and successfully recognized by the Bayesian blocks method.

Zoom In Zoom Out Reset image size

After evaluating the ${F}_{{\rm{B}}}$ scans for the different simulated combinations of background and injected flares, the value of ${F}_{{\rm{B}}}$ to be used in this analysis was set to 5. Using this value, the rate of fake flares drops significantly, while the rate of finding the injected flares is still high. Figure 10 shows an example of the ${F}_{{\rm{B}}}$ scan, and in Figure 11 an example of the denoising method applied to a real light curve for one of the candidate sources with the chosen value of ${F}_{{\rm{B}}}=5$ is shown.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At this point some selection criteria were applied on the Fermi-LAT monitored list of sources, with the aim of selecting light curves with flaring behavior and significant enhancement of the flux over the average level. The first criterion was that the denoised light-curve flux must reach above ${10}^{-6}\;\mathrm{photons}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ during IC-86I.⁵⁸ The second selection criterion aims at identifying denoised light curves exhibiting significant variations in time. In order to quantify this, an 11 day running mean⁵⁹ was calculated for each of the candidate sources. The maximum spread of the running mean was then divided by the mean over the entire 3 yr data-taking period, resulting in a measure of the maximum relative time variation. Only sources for which this maximum relative time variation was greater than 0.5 were selected.

Table 2.  Location and Best-fit Parameters (Number of Signal Events, Spectral Index, Mean Time, and Width Parameter of Gaussian Fit) of the Most Significant Point (Smallest Pre-trial p value) Found in the Data Sample and the Post-trial p value (i.e., Final Significance)

Data Set	R.A.	Decl.	${\hat{n}}_{{\rm{s}}}$	$\hat{\gamma }$	${\hat{T}}_{0}$ (MJD)	${\hat{\sigma }}_{T}$ (days)	p value Pre-trial	p value Post-trial
IC-59	2135	−025	14.5	3.9	55,259	5.5	$2.04\times {10}^{-7}$	1.4%
IC-79	34345	−3165	7.2	3.95	55,466	1.8	$1.07\times {10}^{-5}$	66.0%
IC-86I	23595	4295	13.1	2.0	55,882	7.57	$1.06\times {10}^{-5}$	63.0%

Download table as:  ASCIITypeset image

After selecting a set of potentially interesting neutrino sources, the statistical approach described in Section 4 was applied. The general form of the likelihood function given by Equation (1) was used, and the signal pdf as defined in Equation (2) was obtained using the denoised Fermi-LAT light curves. For each candidate source the likelihood function was maximized with respect to the number of signal events ${n}_{{\rm{s}}}$, the power-law index γ, the time lag ${D}_{{\rm{t}}}$ and flux threshold ${f}_{\mathrm{th}}$. The time lag parameter allows for a time offset between the photon light curve and the neutrino pdf of up to ±0.5 days. Since neutrinos are expected to be produced in the high-activity states of the sources, the flux threshold was varied during the maximization procedure. Once the threshold changes, the time pdf is redefined, setting it equal to zero below the threshold and normalizing to unity what was left above the threshold. The test statistic for this search is given by the maximum likelihood ratio:

$$\begin{eqnarray}&&{TS}=-2\mathrm{log}\left[\displaystyle \frac{{\mathcal{L}}({n}_{{\rm{s}}}=0)}{{\mathcal{L}}({\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{D}}_{{\rm{t}}},{\hat{f}}_{\mathrm{th}})}\right],\end{eqnarray} \tag{ 9 }$$

where ${\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{D}}_{{\rm{t}}},{\hat{f}}_{\mathrm{th}}$ are the best-fit values for the number of signal events, the power-law index, the time lag, and time pdf threshold.

6.2. Results

Using the selection criteria above, the list of sources in Table 3 was selected. The most significant deviation from the background-only hypothesis was observed for the quasar PKS 2142–75 at (R.A., decl.) = (3268, −756). To evaluate the post-trial p value, time-scrambled samples of the IceCube events were generated, and the analysis was repeated on them for all the selected sources. Then the most significant result for each of the scrambled data sets was identified and its significance compared to the p value for our most significant IceCube result. In 77% of the scrambled sets the p value was equal to or smaller than the most significant p value observed in the non-scrambled data, and therefore it is well compatible with background fluctuations.

Table 3.  Results of the Triggered Multi-messenger Flare Search

Name	Decl. (°)	R.A. (°)	p value	${\hat{n}}_{{\rm{s}}}$	${\hat{\gamma }}_{{\rm{s}}}$	${\hat{D}}_{{\rm{t}}}$ (days)	${\hat{f}}_{\mathrm{th}}$ ($\mathrm{ph}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$)	Duration (days)	$\mathrm{log}\left(\displaystyle \frac{{E}_{\mathrm{min}}}{\mathrm{GeV}}\right)$	$\mathrm{log}\left(\displaystyle \frac{{E}_{\mathrm{max}}}{\mathrm{GeV}}\right)$	Fluence Limit (GeV cm−2)
PKS 2142–75	−75.6	326.8	0.02	1.9	3.95	−0.40	1.32e-06	5	5.0	7.8	8.10
PKS 0235–618	−61.6	39.2	...	0	...	...	...	...	...	...	...
PMN J1038–5311	−53.2	159.7	0.47	0.97	2.25	0.50	0	1075	5.2	7.8	9.13
Fermi J1717–5156	−51.9	259.4	0.50	0.32	3.95	0.15	0	1075	5.2	7.8	8.38
PKS 2326–502	−49.9	352.3	...	0	...	...	...	...	...	...	...
PKS 1424–41	−42.1	217.0	0.34	3.02	2.95	−0.37	8.15e-07	910	5.2	7.9	9.73
PKS 0920–39	−40.0	140.7	0.10	6.38	3.25	−0.44	7.93e-07	684	5.1	7.8	11.60
PKS 0426–380	−37.9	67.2	...	0	...	...	...	...	...	...	...
PKS 0402–362	−36.1	61.0	...	0	...	...	...	...	...	...	...
PMN J2250–2806	−28.1	342.7	...	0	...	...	...	...	...	...	...
PKS 2255–282	−28.0	344.5	...	0	...	...	...	...	...	...	...
PKS 1244–255	−25.8	191.7	0.08	0.98	3.22	0.03	1.73e-06	2	5.1	7.9	3.51
PKS 1622–253	−25.5	246.4	0.49	0.54	2.75	−0.20	2.28e-06	82	5.1	7.9	2.37
PMN J1626–2426	−24.4	246.8	0.30	3.28	2.15	−0.08	1.98e-06	269	5.2	7.8	3.58
PKS 0454–234	−23.4	74.3	0.16	1.40	2.05	−0.30	2.38e-06	1	5.1	7.8	2.92
PKS 1830–211	−21.1	278.4	...	0	...	...	...	...	...	...	...
PMN J2345–1555	−15.9	356.3	...	0	...	...	...	...	...	...	...
Fermi J1532–1321	−13.4	233.2	0.45	2.56	2.85	0.40	0	1075	4.9	7.8	2.57
PKS 1730–130	−13.1	263.3	...	0	...	...	...	...	...	...	...
PKS 0727–11	−11.7	112.6	0.09	7.74	3.55	−0.06	8.64e-07	1069	4.8	7.8	3.43
PKS 1346–112	−11.5	207.4	0.18	1.37	2.21	−0.50	1.50e-06	2	4.8	7.8	1.49
PKS 1510–089	−8.8	228.2	0.45	1.68	3.56	−0.17	0	1075	4.5	7.7	1.63
3C 279	−5.8	194.0	...	0	...	...	...	...	...	...	...
PKS 2320–035	−3.3	350.9	0.12	2.67	2.50	0.28	1.34e-06	4	3.5	7.3	0.40
PMN J0948+0022	0.4	147.2	0.27	7.12	3.95	0.50	0	1075	3.5	7.2	0.72
3C 273	2.1	187.3	0.26	1.65	1.85	−0.36	1.12e-06	201	3.4	7.0	0.38
PMN J0505+0416	4.3	76.4	...	0	...	...	...	...	...	...	...
J123939+044409	4.7	189.9	...	0	...	...	...	...	...	...	...
OG 050	7.5	83.2	...	0	...	...	...	...	...	...	...
CTA 102	11.7	338.2	0.26	3.58	3.95	0.10	7.14e-07	128	3.2	6.4	0.41
3C 454.3	16.1	343.5	0.28	1.94	3.65	−0.45	1.99e-06	10	3.2	6.3	0.34
OX 169	17.7	325.9	0.20	2.70	3.95	−0.50	9.45e-07	29	3.2	6.2	0.44
OJ 287	20.1	133.7	...	0	...	...	...	...	...	...	...
PKS B1222+216	21.4	186.2	...	0	...	...	...	...	...	...	...
Crab Pulsar	22.0	83.6	...	0	...	...	...	...	...	...	...
4C 28.07	28.8	39.5	...	0	...	...	...	...	...	...	...
Ton 599	29.2	179.9	...	0	...	...	...	...	...	...	...
B2 1520+31	31.7	230.5	0.42	3.75	2.15	0.44	0	1075	3.0	5.8	0.78
B2 1846+32	32.3	282.1	0.37	1.57	3.95	0.37	8.58e-07	81	3.0	5.8	0.43
B2 0619+33	33.4	95.7	...	0	...	...	...	...	...	...	...
1H 0323+342	34.2	51.2	0.47	1.52	2.15	−0.25	5.40e-07	1051	3.0	5.8	0.71
4C 38.41	38.1	248.8	0.23	11.55	2.55	0.50	0	1075	3.0	5.7	1.19
3C 345	39.7	250.4	0.09	4.0	2.18	0.16	1.07e-06	61	3.0	5.7	0.83
NGC 1275	41.5	50.0	0.15	2.67	3.95	−0.40	1.25e-06	24	3.0	5.7	0.60
BL Lac	42.3	330.7	...	0	...	...	...	...	...	...	...
B3 1343+451	44.9	206.4	0.30	0.82	1.55	0.19	1.29e-06	14	3.0	5.6	0.50
4C 49.22	49.5	178.4	0.11	4.72	3.95	0.25	4.29e-07	99	3.0	5.6	0.66
NRAO 676	50.8	330.4	...	0	...	...	...	...	...	...	...
BZU J0742+5444	54.7	115.7	...	0	...	...	...	...	...	...	...
S4 1849+67	67.1	282.3	0.25	11.70	3.75	0.18	2.90e-07	1046	2.9	5.2	1.27
S5 0836+71	70.9	130.4	0.23	3.13	3.95	−0.32	5.96e-07	109	2.9	5.2	0.75
PKS 0716+714	71.4	110.4	...	0	...	...	...	...	...	...	...
S5 1803+78	78.5	270.2	0.19	3.98	3.95	0.12	7.92e-07	26	2.9	5.2	0.95

Note. Besides the names and the equatorial coordinates of the sources, the best-fit values for the likelihood function parameters ${\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{D}}_{{\rm{t}}},{\hat{f}}_{\mathrm{th}}$ are given. The p values are pre-trial. The columns duration (i.e., the total amount of time for which the Fermi-LAT light curve is above the threshold ${\hat{f}}_{\mathrm{th}}$), $\mathrm{log}({E}_{\mathrm{min}})$, and $\mathrm{log}({E}_{\mathrm{max}})$ list the values used for calculating the 90% C.L. fluence upper limit. The fluence upper limits were calculated with time dependence corresponding to the fitted signal in the likelihood (i.e., non-zero only when the Fermi-LAT light curve is above the threshold ${\hat{f}}_{\mathrm{th}}$) and are hence provided only for sources that were positive fluctuations. These limits were calculated for a flux with spectral index of 2 irrespective of the actual best-fit value. For sources with under-fluctuation, the number of fitted signal events ${\hat{n}}_{{\rm{s}}}$ is zero.

Download table as:  ASCIITypeset image

Figure 12 shows, for PKS 2142–75, the best-fit flux threshold, together with the denoised light curve and the IceCube event weights w_i defined in Equation (5). In this figure one can see that the fit prefers a high flux threshold value, therefore reducing the time pdf to be non-zero only in a narrow time interval, leading to a low best-fit signal strength of ${\hat{n}}_{{\rm{s}}}=1.9$

Zoom In Zoom Out Reset image size

In the searches presented in the previous section neutrino emission was assumed to follow the γ-ray photon flux of the source, using light curves extracted from Fermi-LAT data in the energy range from 100 MeV to 300 GeV as templates. A complementary search was performed to cover a potentially interesting group of sources that did not pass the selection criteria in Section 6.1. These are sources that exhibited flares in the TeV range, but which did not show significant activity in the lower energy range covered by the Fermi-LAT light curves. For these flares Astronomer's Telegrams (ATel) were issued by imaging air-Cerenkov telescopes such as H.E.S.S., MAGIC, or VERITAS. As explained in Section 2, such orphan flares, exhibiting TeV emission while lacking emissions in the lower-energy region, may be interesting for hadronic models, if the lack of emission is real and not due to limited exposure of experimental observations (Halzen & Hooper 2005; Reimer et al. 2005; Sahu et al. 2013). This search was performed for the IC-79 and IC-86I periods only since the data from previous periods have already been analyzed (Abbasi et al. 2012a).

7.1. Method

7.2. Results

In Table 4 the candidate sources are listed. These were selected from the reports found in the ATel indicated in the table. Two sources, 1ES 0806+524 during IC-79 and PG 1553+113 during IC-86I, showed a positive fluctuation over the background, but the p values are compatible with the background.

Table 4.  Source Candidates Selected for the "Search for Triggered Flares with Sporadic Coverage"

Season	Source	ATel num.	Time Range in MJD	${\hat{n}}_{{\rm{s}}}$	${\hat{\gamma }}_{{\rm{s}}}$	p value
IC-79	1ES 0806+524	3192	55615 < T < 55618	0.78	3.95	0.24
	HESS J0632+057	3153, 3161	55598 < T < 55602	0	...	...
	1ES 1215+303	3100	55562 < T < 55565	0	...	...
IC-86I	PG 1553+113	4069	56036 < T < 56039	0.8	3.95	0.23
	BL Lacertae	3459	55739 < T < 55742	0	...	...

Note. The best-fit values for the number of signal events ${\hat{n}}_{{\rm{s}}}$ and the spectral index ${\hat{\gamma }}_{{\rm{s}}}$ are listed. The p values are pre-trial; the ellipsis points indicate that the source best fit was for zero signal events.

Download table as:  ASCIITypeset image

X-ray binaries that have also been observed to emit TeV γ-rays are potential candidates for neutrino emission. While evidence from multiwavelength observations mostly favors leptonic emission (Albert et al. 2009; Hinton et al. 2009), the possibility that protons or nuclei are being accelerated in the jets of these binary systems cannot be ruled out. The observation of neutrino emission from these sources would provide clear evidence for the presence of a hadronic component.

Neutrinos could be produced in these binary systems during interactions of accelerated protons in relativistic jets with the atmosphere of the binary star (Reynoso et al. 2008). These jets are narrow and precess with the same time period as the binary system. The neutrino flux at Earth from these sources depends on the orientation of these jets with respect to the atmosphere of the massive star and our line of sight and hence is expected to be high only during a narrow time window of the orbit.

The list of sources considered and the motivations are explained in detail in references (Abbasi et al. 2012c). Three new sources were added to this list, for a total of 10 sources (see Table 5). In the northern sky a recently reported binary, HESS J0632+057 (Hinton et al. 2009), was added, which is a variable point-like source of VHE ($\gt 100$ GeV) γ-rays located in the Galactic plane and is positionally coincident with a Be star. It also emits variably in the radio and X-ray domains and has been found to have a hard X-ray spectrum (Hinton et al. 2009). The periodicity of the X-ray emission is ${\rm{\Omega }}=320\pm 5$ days (Bongiorno et al. 2011). Bearing a close resemblance to the source LS I+61 303, this source has now been confirmed to be a γ-ray binary (Aleksic et al. 2012a). Motivated by the increased sensitivity of IceCube to neutrino sources in the southern sky as a consequence of background rejection techniques introduced recently with the first year of data from the completed IC-86I configuration of the detector (Aartsen et al. 2013b), two sources in the southern sky were also added. LS 5039 is a high-mass X-ray binary that was also the first microquasar to be established as a high-energy γ-ray source (Aharonian et al. 2005). It has been confirmed to have a period of 3.906 ± 0.002 days (Aharonian et al. 2005). The second new source in the southern sky is HESS J1018–589, a gamma-ray binary in the Carina arm region of the galaxy. This source position is coincident with 1FGL J1018.6 reported by the Fermi-LAT Collaboration (Ackermann et al. 2012) and is a source of ($\gt 100$ GeV) γ-rays. Its periodicity has been reported to be 16.58 ± 0.02 days (Abramowski et al. 2012).

Table 5.  Candidate Sources for the Periodic Neutrino Emission Search

Source (Reference)	Period (days)	p value	Flare Duration	Phase	Time-dependent 90% C.L. Upper Limit (GeV−1 cm−2 s−1)	Time-integrated 90% C.L. Upper Limit (GeV−1 cm−2 s−1)
Cygnus X-1 (Brocksopp et al. 1999)	5.599829 ± 0.000016	0.45	0.016	0.81	$2.31\times {10}^{-10}$	$2.33\times {10}^{-9}$
Cygnus X-3 (Abdo et al. 2009)	0.199679 ± 0.000003	0.34	0.05	0.080	$5.05\times {10}^{-10}$	$1.70\times {10}^{-9}$
GRO J0422+32 (Chevalier & Ilovaisky 1996)	0.212140 ± 0.000003	...	...	...	...	$1.78\times {10}^{-9}$
GRS 1915+105 (Neil et al. 2006)	30.8 ± 0.2	0.31	0.17	0.28	$3.33\times {10}^{-10}$	$1.18\times {10}^{-9}$
LSI + 61 303 (Albert et al. 2009)	26.496 ± 0.0028	...	...	...	...	$1.95\times {10}^{-9}$
SS 433 (Kubota et al. 2010)	13.08227 ± 0.00008	...	...	...	...	$6.5\times {10}^{-10}$
XTE J1118+480 (Torres et al. 2004)	0.1699339 ± 0.0000002	...	...	...	...	$1.21\times {10}^{-9}$
HESS J0632+057 (Bongiorno et al. 2011)	320 ± 5	0.087	0.0127	0.70	$4.82\times {10}^{-10}$	$1.37\times {10}^{-9}$
LS 5039 (Aharonian et al. 2005)	3.906 ± 0.002	...	...	...	...	$5.24\times {10}^{-9}$
HESS J1018–589 (Abramowski et al. 2012)	16.58 ± 0.02	...	...	...	...	$9.21\times {10}^{-9}$

Note. The p values are pre-trial. The flare durations given are the widths of the the fitted Gaussian flares, as fractions of the total periods. The time-dependent upper limits are the normalization for an ${E}^{-2}$ power-law flux with time dependence corresponding to the fitted signal in likelihood and are hence provided only for sources that were positive fluctuations. Time-integrated upper limits (Aartsen et al. 2013b) are a factor of 3–8 times lower than reported by analyses performed on data from the IC-22 and IC-40 configurations of the detector (Abbasi et al. 2012c), and the time-dependent upper limits have improved similarly, where comparable.

Download table as:  ASCIITypeset image

8.1. Method

The analysis method used here was previously applied to data from the IC-22 and IC-40 configurations of the IceCube detector and is explained in detail in Abbasi et al. (2012c). It is adapted from the method used in the All-Sky Time Scan (Section 5). Each event within the sample is assigned an event phase ${\phi }_{i}$ based on the arrival time of the event and the period of the source, known mainly from optical measurements. Consequently, signal events coming from different periods of a periodic signal are assigned similar values of ${\phi }_{i}$.

The time dependence of the signal pdf ${{\mathcal{T}}}_{i}^{\mathrm{signal}}$ in Equation (2) is then replaced with a phase dependence ${{\mathcal{P}}}_{i}^{\mathrm{signal}}$ of the form

$$\begin{eqnarray}&&{{\mathcal{P}}}_{i}^{\mathrm{signal}}=\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{\Phi }}}}\mathrm{exp}\left(-\displaystyle \frac{{({\phi }_{i}-{\rm{\Phi }})}^{2}}{2{\sigma }_{{\rm{\Phi }}}^{2}}\right),\end{eqnarray} \tag{ 10 }$$

where ${\phi }_{i}$ is the arrival phase of the ith event and Φ and ${\sigma }_{{\rm{\Phi }}}$ are the position of the peak signal in phase and the width, respectively. These two parameters are fitted to maximize the likelihood. Two more free parameters were part of the fit, the number of signal events ${n}_{{\rm{s}}}$ and the spectral index γ. Thus, the method searches for a statistically significant clustering of high-energy events not only in space but also in phase.

In Abbasi et al. (2012c) this search was constrained to look for a flare of duration larger than $0.02\times {\rm{\Omega }}$, where Ω is the period. This was done in order to avoid the fact that the minimizer tends to prefer shorter flares given a limited time window. The same marginalization term described previously for the All-Sky Time Scan was used in this search to prevent very short flares from dominating the significance, while still allowing a search for short flares that are of physical interest. The test statistic for this search including this modification is thus

$$\begin{eqnarray}&&{TS}=-2\mathrm{log}\left[\displaystyle \frac{1}{\sqrt{2\pi }{\hat{\sigma }}_{{\rm{\Phi }}}}\times \displaystyle \frac{{\mathcal{L}}({n}_{{\rm{s}}}=0)}{{\mathcal{L}}({\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{\sigma }}_{{\rm{\Phi }}},\hat{{\rm{\Phi }}})}\right],\end{eqnarray} \tag{ 11 }$$

where ${\hat{n}}_{{\rm{s}}},{\hat{\gamma }}_{{\rm{s}}},{\hat{\sigma }}_{{\rm{\Phi }}},\hat{{\rm{\Phi }}}$ are the best-fit values and $\frac{1}{\sqrt{2\pi }{\hat{\sigma }}_{{\rm{\Phi }}}}$ is the marginalization term.

The post-trial p value of the most significant observation was estimated by calculating the fraction of time-scrambled data sets in which the most significant fluctuation observed among the 10 sources was more significant than that which was observed in non-scrambled data. As can be seen in Figure 13, the search is more sensitive to flares of very short duration measured as a fraction of the period of the system.

Zoom In Zoom Out Reset image size

As the considered periodic sources are not expected to change their long-term behavior, the performance of this search improves as more data are included. For the source GRS 1915+105, which has large relative uncertainty on the measured period, converting the event times into phases using the reported central value may lead to a smearing out of the actual flare if it happened at a different time within the error. Since the signal is modeled as a Gaussian flare in phase (Abbasi et al. 2012c) and is more sensitive to narrower flares (Figure 13), this could negatively impact the sensitivities and the discovery potentials. The effect is proportional to the relative error on the period times the number of periods within the livetime and is hence most severe for GRS 1915+105. The impact of the period uncertainty can be estimated by recalculating the discovery potential assuming that the true source period differs from the central by the reported uncertainty. In Figure 14 the result of this test is shown, demonstrating that this effect is not sufficient to negate the comparative advantage this search has over time-integrated searches.

Zoom In Zoom Out Reset image size

8.2. Results

The results of the periodic analysis for each of the selected sources are given in Table 5. The most significant observation was from the source HESS J0632+057 with a pre-trial probability of 8.67%. This Gaussian-fitted flare was observed at a phase of 0.702 with a width of ${\sigma }_{{\rm{\Phi }}}=0.012$ in terms of the fraction of the period. Figure 15 shows the fitted Gaussian and the IceCube event weights w_i defined in Equation (5). Cyg X-1, Cyg X-3, and GRS 1915+105 were observed to have flares of probability 0.45, 0.34, and 0.32, respectively. All other sources produced under-fluctuations, indicating that the number of events in the direction of the source was less than or equal to the number expected from background only. The post-trial probability of the fluctuation from HESS J0632+057 was found to be 44.3%, making the observation compatible with the background-only hypothesis.

Zoom In Zoom Out Reset image size

Searches described within this paper have found no evidence for the existence of flaring or periodic neutrino sources. Both the untriggered search looking for neutrino flares anywhere in the sky and the triggered search looking for neutrino flares coinciding with γ-ray flares reported by the Fermi-LAT returned results consistent with the background-only hypothesis. No evidence has been found for periodic neutrino emission by binary systems either. These analyses include data from the first year of operation of the completed IceCube detector, taken between 2011 May and 2012 May. IceCube will continue to run in this configuration for the foreseeable future, and the resultant signal buildup will increase the sensitivity to steady point sources of neutrinos (Aartsen et al. 2013b). A similar improvement is expected in the sensitivity of the detector toward periodic neutrino signals. For flare searches, on the other hand, additional years of data with the full detector improve the chances to see rarer (rather than weaker) neutrino emission from outbursts, which we have not been fortunate enough to witness yet.

In this paper we have demonstrated the viability of long-term monitoring of sources of interest triggered by multiwavelength information from other experiments. As the detector matures, detector operations, data acquisition, and processing will become more automated, and we will soon be able to carry out this monitoring near-real-time—reducing the time delay between the trigger and the results.

We acknowledge the support from the following agencies: U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation—Physics Division, University of Wisconsin Alumni Research Foundation, the Grid Laboratory Of Wisconsin (GLOW) grid infrastructure at the University of Wisconsin−Madison, the Open Science Grid (OSG) grid infrastructure; U.S. Department of Energy, and National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; Natural Sciences and Engineering Research Council of Canada, WestGrid and Compute/Calcul Canada; Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation, Sweden; German Ministry for Education and Research (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Research Department of Plasmas with Complex Interactions (Bochum), Germany; Fund for Scientific Research (FNRS-FWO), FWO Odysseus programme, Flanders Institute to encourage scientific and technological research in industry (IWT), Belgian Federal Science Policy Office (Belspo); University of Oxford, United Kingdom; Marsden Fund, New Zealand; Australian Research Council; Japan Society for Promotion of Science (JSPS); the Swiss National Science Foundation (SNSF), Switzerland; National Research Foundation of Korea (NRF); Danish National Research Foundation, Denmark (DNRF).