Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert

See allHide authors and affiliations

Science  12 Jul 2018:
DOI: 10.1126/science.aat2890


A high-energy neutrino event detected by IceCube on 22 September 2017 was coincident in direction and time with a gamma-ray flare from the blazar TXS 0506+056. Prompted by this association, we investigated 9.5 years of IceCube neutrino observations to search for excess emission at the position of the blazar. We found an excess of high-energy neutrino events, with respect to atmospheric backgrounds, at that position between September 2014 and March 2015. Allowing for time-variable flux, this constitutes 3.5σ evidence for neutrino emission from the direction of TXS 0506+056, independent of and prior to the 2017 flaring episode. This suggests that blazars are identifiable sources of the high-energy astrophysical neutrino flux.

The origin of the highest-energy cosmic rays is believed to be extragalactic (1), but their acceleration sites remain unidentified. High-energy neutrinos are expected to be produced in or near the acceleration sites when cosmic rays interact with matter and ambient light, producing charged mesons that decay into neutrinos and other particles. Unlike cosmic rays, neutrinos can travel through the Universe unimpeded by interactions with other particles and undeflected by magnetic fields, providing a means to identify and study the extreme environments producing cosmic rays (2). Blazars, a class of active galactic nuclei with powerful relativistic jets pointed close to our line of sight (3), are prominent candidate sources of such high-energy neutrino emission (49). The electromagnetic emission of blazars is observed to be highly variable on time scales from minutes to years (10).

The IceCube Neutrino Observatory (11) is a high-energy neutrino detector occupying an instrumented volume of 1 km3 within the Antarctic ice sheet at the Amundsen-Scott South Pole Station. The detector consists of an array of 86 vertical strings, nominally spaced 125 m apart and descending to a depth of approximately 2450 m in the ice. The bottom 1 km of each string is equipped with 60 optical sensors that record Cherenkov light emitted by relativistic charged particles passing through the optically transparent ice. When high-energy muon neutrinos interact with the ice, they can create relativistic muons that travel many kilometers, creating a track-like series of Cherenkov photons recorded when they pass through the array. This allows the reconstruction of the original neutrino direction with a median angular uncertainty of 0.5° for a neutrino energy of ~30 TeV (or 0.3° at 1 PeV) (12, 13).

IceCube discovered the existence of a diffuse flux of high-energy astrophysical neutrinos in 2013 (14, 15). Measurements of the energy spectrum have since been refined (16, 17), indicating that the neutrino spectrum extends above several PeV. However, analyses of neutrino observations have not succeeded in identifying individual sources of high-energy neutrinos (12, 18). This suggests that the sources are distributed across the sky and that even the brightest individual sources contribute only a small fraction of the total observed flux.

Recently, the detection of a high-energy neutrino by IceCube, together with observations in gamma rays and at other wavelengths, indicates that a blazar, TXS 0506+056, located at right ascension (RA) 77.3582° and declination (Dec) +5.69314° (J2000 equinox) (19) may be an individually identifiable source of high-energy neutrinos (20). The neutrino-candidate event, IceCube-170922A, was detected on 22 September 2017, selected by the Extremely High Energy (EHE) online event filter (21), and reported as a public alert (22). EHE alerts are currently sent at a rate of about four per year, and are based on well-reconstructed, high-energy muon-track events. The selection threshold is set so that approximately half of the events are estimated to be astrophysical neutrinos, the rest being atmospheric background events. After the alert was sent, further studies refined the directional reconstruction, with best-fitting coordinates of RA Embedded Image and Dec Embedded Image (degrees, J2000, 90% containment region). The most probable neutrino energy was estimated to be 290 TeV, with a 90% confidence level lower limit of 183 TeV (20).

It was soon determined that the direction of IceCube-170922A was consistent with the location of TXS 0506+056 and coincident with a state of enhanced gamma-ray activity observed since April 2017 (23) by the Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope (24). Follow-up observations of the blazar led to the detection of gamma rays with energies up to 400 GeV by the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) Telescopes (25, 26). IceCube-170922A and the electromagnetic observations are described in detail in (20). The significance of the spatial and temporal coincidence of the high-energy neutrino and the blazar flare is estimated to be at the 3σ level (20). On the basis of this result, we consider the hypothesis that the blazar TXS 0506+056 has been a source of high-energy neutrinos beyond that single event.

Searching for neutrino emission

IceCube monitors the whole sky and has maintained essentially continuous observations since 5 April 2008. Searches for neutrino point sources using two model-independent methods, a time-integrated and a time-dependent unbinned maximum likelihood analysis, have previously been published for the data collected between 2008 and 2015 (12, 18, 27). Here, we analyze the same 7-year data sample supplemented with additional data collected from May 2015 until October 2017 (21). The data span 9.5 years and consist of six distinct periods, corresponding to changing detector configurations, data-taking conditions, and improved event selections (Table 1).

Table 1 IceCube neutrino data samples.

Six data-taking periods make up the full 9.5-year data sample. Sample numbers correspond to the number of detector strings that were operational. During the first three periods, the detector was still under construction. The last three periods correspond to different data-taking conditions and/or event selections with the full 86-string detector.

View this table:

The northern sky, where TXS 0506+056 is located, is observed through Earth by IceCube. Approximately 70,000 neutrino-induced muon tracks are recorded each year from this hemisphere of the sky after passing the final event selection criteria. Fewer than 1% of these events originate from astrophysical neutrinos; the vast majority are background events caused by neutrinos of median energy ~1 TeV created in cosmic ray interactions in the atmosphere over other locations on Earth. However, for an astrophysical muon-neutrino flux where the differential number of neutrinos with energy E scales as dN/dE ~ E–2, the distribution of muon energies is different than for the background atmospheric neutrino flux, which scales as ~E–3.7 (17). This allows for further discriminating power in point source searches besides directional-only excesses.

A high-significance point source detection (12, 18) can require as few as two or three, or as many as 30, signal events to stand out from the background, depending on the energy spectrum and the clustering of events in time. To search for a neutrino signal at the coordinates of TXS 0506+056, we apply the standard time-integrated analysis (28) and time-dependent analysis (29) that have been used in past searches (12, 18, 27). The time-integrated analysis uses an unbinned maximum likelihood ratio method to search for an excess number of events consistent with a point source at a specified location, given the angular distance and angular uncertainty of each event. Energy information is included in the definition of the likelihood, assuming a power-law energy spectrum E–γ , with the spectral index γ as a fitted parameter. The model parameters are correlated and are expressed as a pair, (Φ100, γ), where Φ100 is the flux normalization at 100 TeV. The time-dependent analysis uses the same formulation of the likelihood but searches for clustering in time as well as space by introducing an additional time profile. It is performed separately for two different generic profile shapes: a Gaussian-shaped time window and a box-shaped time window. Each analysis varies the central time of the window, T0, and the duration TW (from seconds to years) of the potential signal to find the four parameters (Φ100, γ, T0, TW) that maximize the likelihood ratio, which is defined as the test statistic TS. (For the Gaussian time window, TW represents twice the standard deviation.) The test statistic includes a factor that corrects for the look-elsewhere effect arising from all of the possible time windows that could be chosen (30).

For each analysis method (time-integrated and time-dependent), a robust significance estimate is obtained by performing the identical analysis on trials with randomized datasets. These are produced by randomizing the event times and recalculating the RA coordinates within each data-taking period. The resultant P value is defined as the fraction of randomized trials yielding a value of TS greater than or equal to the one obtained for the actual data.

Because the detector configuration and event selections changed as shown in Table 1, the time-dependent analysis is performed by operating on each data-taking period separately. (A flare that spans a boundary between two periods could be partially detected in either period, but with reduced significance.) An additional look-elsewhere correction then needs to be applied for a result in an individual data segment, given by the ratio of the total 9.5-year observation time to the observation time of that data segment (30).

Neutrinos from the direction of TXS 0506+056

The results of the time-dependent analysis performed at the coordinates of TXS 0506+056 are shown in Fig. 1 for each of the six data periods. One of the data periods, IC86b from 2012 to 2015, contains a significant excess, which is identified by both time-window shapes. The excess consists of 13 ± 5 events above the expectation from the atmospheric background. The significance depends on the energies of the events, their proximity to the coordinates of TXS 0506+056, and their clustering in time. This is illustrated in Fig. 2, which shows the time-independent weight of individual events in the likelihood analysis during the IC86b data period.

Fig. 1 Time-dependent analysis results.

The orange curve corresponds to the analysis using the Gaussian-shaped time profile. The central time T0 and width TW are plotted for the most significant excess found in each period, with the P value of that result indicated by the height of the peak. The blue curve corresponds to the analysis using the box-shaped time profile. The curve traces the outer edge of the superposition of the best-fitting time windows (durations TW) over all times T0, with the height indicating the significance of that window. In each period, the most significant time window forms a plateau, shaded in blue. The large blue band centered near 2015 represents the best-fitting 158-day time window found using the box-shaped time profile. The vertical dotted line in IC86c indicates the time of the IceCube-170922A event.

Fig. 2 Time-independent weight of individual events during the IC86b period.

Each vertical line represents an event observed at the time indicated by calendar year (top) or MJD (bottom). Overlapping lines are shifted by 1 to 2 days for visibility. The height of each line indicates the event weight: the product of the event’s spatial term and energy term in the unbinned likelihood analysis evaluated at the location of TXS 0506+056 and assuming the best-fitting spectral index γ = 2.1 (30). The color for each event indicates an approximate value in units of TeV of the reconstructed muon energy (muon energy proxy), which the analysis compares with expected muon energy distributions under different hypotheses. [A distribution for the true neutrino energy of a single event can also be inferred from the event’s muon energy (30).] The dashed curve and the solid bracket indicate the best-fitting Gaussian and box-shaped time windows, respectively. The distribution of event weights and times outside of the best-fitting time windows is compatible with background.

The Gaussian time window is centered at 13 December 2014 [modified Julian day (MJD) 57004] with an uncertainty of ±21 days and a duration TW = Embedded Image days. The best-fitting parameters for the fluence J100 = ∫Φ100(t)dt and the spectral index are given by E2J100 = Embedded Image TeV cm–2 at 100 TeV and γ = 2.1 ± 0.2, respectively. The joint uncertainty on these parameters is shown in Fig. 3 along with a skymap showing the result of the time-dependent analysis performed at the location of TXS 0506+056 and in its vicinity during the IC86b data period.

Fig. 3 Time-dependent analysis results for the IC86b data period (2012–2015).

(A) Change in test statistic, ΔTS, as a function of the spectral index parameter γ and the fluence at 100 TeV given by E2J100. The analysis is performed at the coordinates of TXS 0506+056, using the Gaussian-shaped time window and holding the time parameters fixed (T0 = 13 December 2014, TW = 110 days). The white dot indicates the best-fitting values. The contours at 68% and 95% confidence level assuming Wilks’ theorem (36) are shown in order to indicate the statistical uncertainty on the parameter estimates. Systematic uncertainties are not included. (B) Skymap showing the P value of the time-dependent analysis performed at the coordinates of TXS 0506+056 (cross) and at surrounding locations. The analysis is performed on the IC86b data period, using the Gaussian-shaped time window. At each point, the full fit for (Φ, γ, T0, TW) is performed. The P value shown does not include the look-elsewhere effect related to other data periods. An excess of events is detected, consistent with the position of TXS 0506+056.

The box-shaped time window is centered 13 days later with duration TW = 158 days (from MJD 56937.81 to MJD 57096.21, inclusive of contributing events at boundary times). For the box-shaped time window, the uncertainties are discontinuous and not well defined, but the uncertainties for the Gaussian window show that it is consistent with the box-shaped time window fit. Despite the different window shapes, which lead to different weightings of the events as a function of time, both windows identify the same time interval as significant. For the box-shaped time window, the best-fitting parameters are similar to those of the Gaussian window, with fluence at 100 TeV and spectral index given by E2J100 = Embedded Image TeV cm–2 and γ = 2.2 ± 0.2. This fluence corresponds to an average flux over 158 days of Φ100 = Embedded Image TeV–1 cm–2 s–1.

When we estimate the significance of the time-dependent result by performing the analysis at the coordinates of TXS 0506+056 on randomized datasets, we allow in each trial a new fit for all the parameters: Φ100, γ, T0, TW. We find that the fraction of randomized trials that result in a more significant excess than the real data is 7 × 10–5 for the box-shaped time window and 3 × 10–5 for the Gaussian time window. This fraction, once corrected for the ratio of the total observation time to the IC86b observation time (9.5 years/3 years), results in P values of 2 × 10–4 and 10–4, respectively, corresponding to 3.5σ and 3.7σ. Because there is no a priori reason to prefer one of the generic time windows over the other, we take the more significant one and include a trial factor of 2 for the final significance, which is then 3.5σ.

Outside the 2012–2015 time period, the next most significant excess is found using the Gaussian window in 2017 and includes the IceCube-170922A event. This time window is centered at 22 September 2017 with duration TW = 19 days, γ = 1.7 ± 0.6, and fluence E2J100 = Embedded Image TeV cm–2 at 100 TeV. No other event besides the IceCube-170922A event contributes significantly to the best fit. As a consequence, the uncertainty on the best-fitting window location and width spans the entire IC86c period, because any window containing IceCube-170922A yields a similar value of the test statistic. Following the trial correction procedure for different observation periods as described above, the significance of this excess is 1.4σ. If the IceCube-170922A event is removed, no excess remains during this time period. This agrees with the result of the rapid-response analysis (31) that is part of the IceCube alert program, which found no other potential astrophysical neutrinos from the same region of the sky during ±7 days centered on the time of IceCube-170922A.

We performed a time-integrated analysis at the coordinates of TXS 0506+056 using the full 9.5-year data sample. The best-fitting parameters for the flux normalization and the spectral index are Φ100 = Embedded Image TeV–1 cm–2 s–1 and γ = 2.0 ± 0.3, respectively. The joint uncertainty on these parameters is shown in Fig. 4A. The P value, based on repeating the analysis at the same coordinates with randomized datasets, is 0.002% (4.1σ), but this is an a posteriori significance estimate because it includes the IceCube-170922A event, which motivated performing the analysis at the coordinates of TXS 0506+056. An unbiased significance estimate including the event would need to take into account the look-elsewhere effect related to all other possible directions in the sky that could be analyzed. It is expected that there will be two or three directions somewhere in the northern sky with this significance or greater, resulting from the chance alignment of neutrinos (12). Here, we are interested in determining whether there is evidence of time-integrated neutrino emission from TXS 0506+056 besides the IceCube-170922A event.

Fig. 4 Time-integrated analysis results.

As in Fig. 3A, but for the time-integrated analysis of TXS 0506+056 using (A) the full 9.5-year sample (2008–2017), and (B) the 7-year sample (2008–2015).

If we remove the final data period IC86c, which contains the event, and perform the analysis again using only the first 7 years of data, we find best-fitting parameters that are nearly unchanged: Φ100 = Embedded Image TeV–1 cm–2 s–1 and γ = 2.1 ± 0.3, respectively. The joint uncertainty on these parameters is shown in Fig. 4B. The P value, using only the first 7 years of data, is 1.6% (2.1σ), based on repeating the analysis at the same coordinates with randomized datasets. These results indicate that the time-integrated fit is dominated by the same excess as found in the time-dependent analysis above, having similar values for the spectral index and total fluence (E2J100 = 2.0 × 10–4 TeV cm–2 at 100 TeV over the 7-year period). This excess is not significant in the time-integrated analysis because of the additional background during the rest of the 7-year period.

Blazars as neutrino sources

The signal identified during the 5-month period in 2014–2015 consists of an estimated 13 ± 5 muon-neutrino events that are present in addition to the expected background. The analysis is unbinned, but the mean background at the declination of TXS 0506+056 is useful for comparison purposes; it is 5.8 events in a search bin of radius 1° during a 158-day time window. (We use the duration of the box-shaped time window result for convenience to calculate averages during the flare.) The significance of the excess is due to both the number of events and their energy distribution, with higher-energy events increasing the significance and leading to the best-fitting spectral index of 2.1, in contrast to the lower-energy atmospheric neutrino background with spectral index ~3.7. At this declination in the sky, the 68% central energy range in which IceCube is most sensitive to point sources with E–2.1 spectra is between 32 TeV and 3.6 PeV. Assuming that the muon-neutrino fluence (E2J100 = Embedded Image TeV cm–2) is one-third of the total neutrino fluence, then the all-flavor neutrino energy fluence is Embedded Image erg cm–2 over this energy range. With the recent measurement (32) of the redshift of TXS 0506+056 as z = 0.3365 ± 0.0010, this energy fluence implies that the isotropic neutrino luminosity is Embedded Image erg s–1 averaged over 158 days. This is higher than the isotropic gamma-ray luminosity during the same period, which is similar to the long-term luminosity between 0.1 GeV and 100 GeV of 0.28 × 1047 erg s–1 averaged over all Fermi-LAT observations of TXS 0506+056 (20). Gamma rays are expected to be produced in the same processes that produce neutrinos—for example, when accelerated protons interact with ambient lower-energy photons near the source, producing both neutral pions (which decay to gamma rays) and charged pions (which decay to neutrinos and leptons). A higher luminosity in neutrinos than in gamma rays could imply that a substantial fraction of the gamma rays related to the neutrino production are either absorbed or arriving at energies above or below the Fermi-LAT energy band.

Although TXS 0506+056 is a bright object in gamma rays, it was not previously singled out as a predicted neutrino source. In the third catalog of active galactic nuclei detected by Fermi-LAT (33) listing 1773 objects (including those at low galactic latitudes), TXS 0506+056 is among the 50 brightest objects, with an average flux between 1 GeV and 100 GeV of 6.5 (±0.2) × 10–9 photons cm–2 s–1. Its measured redshift now makes it one of the most luminous objects known out to the same distance, more than an order of magnitude more luminous than nearby blazars such as Markarian 421, Markarian 501, and 1ES 1959+650. With respect to these objects, an important observational distinction is the favorable declination of TXS 0506+056. As the neutrino-nucleon interaction cross section grows with energy, absorption in Earth becomes considerable for neutrinos above ~100 TeV. IceCube is most sensitive to high-energy neutrinos from sources at declinations near the equatorial plane, which is viewed along the horizon from the South Pole. The blazars mentioned above are at more northern declinations, and the likelihood that a neutrino with energy of ~300 TeV from one of these will be absorbed while traversing Earth is three to five times the likelihood that it will reach the detector. The explanation for why TXS 0506+056 is the first blazar associated with a significant neutrino excess may therefore depend on the combination of its intrinsic properties and the observational properties of the detector.

IceCube recently published (34) a search for neutrino emission from the blazars in the second catalog of active galactic nuclei detected by Fermi-LAT (35), constraining their contribution to the diffuse astrophysical neutrino flux under different model assumptions. An upper limit of 27% was found assuming the diffuse flux that is fit between 10 TeV and 100 TeV with a soft E–2.5 spectrum (16). For an E–2 spectrum compatible with the diffuse flux fit above ~200 TeV (17), the upper limit is between 40% and 80%. The allowed contribution by blazars as a population is larger, because it would include the contribution of fainter and more distant blazars not yet resolved in the catalog. Averaged over 9.5 years, the neutrino flux of TXS 0506+056 by itself corresponds to 1% of the astrophysical diffuse flux and is fully compatible with the previous blazar catalog results.

The evidence presented above supports the hypothesis presented in (20) that the blazar TXS 0506+056 is a high-energy neutrino source. The 3.5σ evidence for neutrino emission during the 5-month period in 2014–2015 is statistically independent of the evidence presented in (20). The analysis of the IceCube-170922A event in (20) relies on correlation of a single neutrino with electromagnetic activity, whereas the analysis presented here relies only on self-correlation of multiple neutrinos. The coincidence of an IceCube alert with a flaring blazar, combined with a neutrino flare from the same object in archival IceCube data, pinpoints a likely source of high-energy cosmic rays.

Supplementary Materials

Materials and Methods

Figs. S1 to S6

IceCube Collaboration Author List

References and Notes

  1. See supplementary materials.
Acknowledgments: Funding: We acknowledge support from the following agencies. United States: National Science Foundation (NSF) Office of Polar Programs, NSF 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), U.S. 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 the astroparticle physics computational facility at Marquette University. Belgium: Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, 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). 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). Author contributions: The IceCube Collaboration designed, constructed, and now operates the IceCube Neutrino Observatory. Data processing and calibration, Monte Carlo simulations of the detector and of theoretical models, and data analyses were performed by a large number of collaboration members, who also discussed and approved the scientific results presented here. The manuscript was reviewed by the entire collaboration before publication, and all authors approved the final version. Competing interests: There are no competing interests to declare. Data and materials availability: Data and resources are available from the IceCube data archive at For each sample listed in Table 1, these include the events within 3° of the TXS 0506+056 source coordinates, neutrino effective areas, background rates, and other supporting information in machine-readable formats.
View Abstract


Navigate This Article