Search for Astronomical Neutrinos from Blazar TXS 0506+056 in Super-Kamiokande

TXS 0506+056 is a BL Lac object type blazar (redshift z = 0.3365 ± 0.0010; Paiano et al. 2018) and located at R.A. = 773582 and decl. = +56931 (J2000 equinox; Massaro et al. 2015). The IceCube Neutrino Observatory (Aartsen et al. 2017) detected a high-energy neutrino event with an estimated energy of 290 TeV on 2017 September 22 at 20:54:30.43 Coordinated Universal Time (IceCube-170922A), the arrival direction of which coincides with the location of TXS 0506+056 (Aartsen et al. 2018a, 2018b). Within a minute of detection, this event's information shared via the Gamma-ray Coordinate Network (GCN)⁵³ and follow-up observations over a wide range of energies were carried out by several observatories. According to the Fermi All-Sky Variability Analysis (Abdollahi et al. 2017), TXS 0506+056 brightened in the GeV band starting in 2017 April (Tanaka et al. 2017). Fermi's Automated Science Processing also found a gamma-ray flare from this source years before. Subsequently, the IceCube collaboration additionally reported a possible neutrino event excess from this blazar in older data between 2014 September and 2015 March (Aartsen et al. 2018b). Coincidence between the neutrino arrival direction and the blazar location as well as timing correlated with the associated gamma-ray flare suggest that the observed neutrinos originated from the blazar and strongly motivate searches for neutrinos in the other energy regions.

Super-Kamiokande (SK; Fukuda et al. 2003) is a large water Cerenkov detector located 1000 m underground (2700 m.w.e.) in the Kamioka-mine, Gifu Prefecture, Japan. It is a cylindrical detector, 39.3 m in diameter and 41.4 m in height and contains 50 kilotonnes of ultra-pure water as neutrino target. The tank is optically separated into an inner detector (ID) and an outer detector (OD) by a structure placed ∼2 m from the tank wall. More than 11,000 twenty-inch photomultiplier tubes (PMTs) in the ID are used to observe the pattern and amount of Cerenkov photons emitted by charged particles produced by neutrino interactions in the water. The OD is primarily used as a veto and has 1885 eight-inch PMTs and is covered with a Tyvek sheet to enhance the light reflection from Cerenkov photons. Super-Kamiokande is primarily sensitive to particle interactions in the energy region of several MeV to a few tens of TeV. In particular, SK observes atmospheric neutrinos above several 10 MeV at a rate of ∼8 events per day in a 22.5 kiloton fiducial volume within the ID (Jiang et al. 2019). The overburden of the mountain above the detector reduces the cosmic-ray muon rate at the detector by a factor of ∼10⁵ compared to that at the surface. Such backgrounds are almost completely eliminated by anticoincidence of the ID and OD, reducing non-neutrino events to less than 1% of the final data sample.

The SK experiment has been operated since 1996 April and has made observations in four distinct phases known as SK-I, SK-II, SK-III, and SK-IV. The first phase, SK-I, lasted from 1996 April to 2001 July (1489.2 livetime days) with 40% photocoverage of the ID using 11,146 PMTs. In 2001 November, half of those PMTs were lost to an accident following detector maintenance, so the SK-II phase was operated from 2002 December to 2005 October (798.6 livetime days) with a reduced photocoverage of 19% (5182 PMTs). After replacing the missing PMTs in 2006 April, the SK-III period operated with the full photocoverage (11,129 PMTs) until 2008 September (518.1 livetime days). New front-end electronics were installed immediately thereafter to start the SK-IV phase. This period lasted until 2018 June, when refurbishment work ahead of a detector upgrade began. Though the detector is currently running as SK-V, in this Letter neutrino data from SK-I to SK-IV through 2018 February corresponding to 5924.4 live days are used for analysis.

The present analysis utilizes the Super-Kamiokande neutrino data with more than 100 MeV of visible energy, divided into three classes depending upon the event topology. In the fully contained (FC) and partially contained (PC) event samples, the neutrino interaction is reconstructed within the ID using Cerenkov rings produced by its daughter particles. An event where all daughter particles stop inside the ID is classified as FC and those where at least one particle exits the ID and deposits energy in the OD is classified as PC. Upward-going muon (UPMU) events are observed when energetic muons produced by muon-neutrino interactions with the rock surrounding the detector penetrate the ID from below its horizon. Since similar downward-going neutrino events suffer from a large amount of cosmic-ray muon backgrounds from above the detector, the event direction of the UPMU sample is restricted to be upward. There is no such restriction on FC and PC events. Both electron neutrinos (ν_e, ${\bar{\nu }}_{e}$) and muon neutrinos (ν_μ, ${\bar{\nu }}_{\mu }$) are observed in the FC sample, while only muon neutrinos populate the PC and UPMU samples.

Events with vertices inside the fiducial volume, defined as the region in the ID more than 2.0 m from any wall, are selected for the FC and PC samples. Separation between the FC and PC samples is determined by the number of effective PMT hits in OD; the FC sample requires fewer than 16 hits (10 hits for SK-I). Events rejected by this cut are classified as PC. UPMU events where the muon passes through the detector (upward through-going muons) as well as events where it stops in the detector (upward stopping muons) are included in this analysis. Upward through-going muon events are required to have a muon track longer than 7.0 m and stopping events are required to have a muon momentum greater than 1.6 GeV. Both criteria ensure that the reconstructed muon direction is from below the horizon. Details of the event selection can be found in Ashie et al. (2005).

In order to estimate the atmospheric neutrino background for this search, a 500 yr equivalent Monte Carlo (MC) simulation of each SK phase has been used. The NEUT interaction generator (Hayato 2009) has been adopted for interactions in water, and a detector simulation based on the Geant3 (Brun et al. 1987) framework is employed for tracking secondary particles and simulating the detector response. Additional simulation details are presented in Abe et al. (2018).

The arrival direction of neutrinos is determined by reconstruction of Cerenkov rings in the ID, and the reconstruction quality typically depends on the number of such rings and the energy of the initial neutrino interaction. In order to use events with sufficient angular resolution for association with the blazar direction in this search, an additional cut on the observed energy is applied, 5.1 GeV for FC events and 1.8 GeV for PC events. This ensures that the angular deviation of the reconstructed direction from the truth is within 10° for more than 68% of these events. Since UPMU events originate from neutrinos with higher energy than other categories, their arrival direction is estimated with higher accuracy. Therefore no additional restriction on the UPMU energy is used as more than 77% of events are reconstructed within 5° of the true arrival direction.

The analysis described below searches for a possible neutrino excess from the blazar by first counting the number of neutrinos in an angular region around the direction to the assumed source. Then the number of events coming from an alternative direction are studied to check the consistency of the observation and background predictions. Finally, a simple statistical method is used to test for local increases in the event rate coming from the blazar direction.

Figure 1 shows a sky map of the reconstructed arrival direction of selected neutrino events around the blazar direction for all samples. Based on the angular resolution of the reconstructed direction in each sample, we defined the search region to be 10° (5°) around the blazar location for FC and PC (UPMU; Abe et al. 2017). There were 18 FC, 29 PC, and 20 UPMU events observed in these regions during SK I–IV. Note that neutrino events observed in the FC sample include both electron-neutrino (ν_e and ${\bar{\nu }}_{e}$) and muon-neutrino (ν_μ and ${\bar{\nu }}_{\mu }$) interactions.

Zoom In Zoom Out Reset image size

In order to quantitatively study possible event excesses above atmospheric neutrino backgrounds, MC is used to predict the event rate in the search region of each sample. The atmospheric neutrino event rate depends on the arrival direction because the thickness of the atmosphere and the neutrino oscillation probability change with zenith angle (in detector coordinates). The thicker the atmosphere, the higher the probability that atmospheric neutrinos will be generated. For downward or horizontally produced neutrinos the path length to the detector is relatively short and the effect of neutrino oscillations is reduced. Consequently, since the zenith angle is related to decl., the event rate also varies with decl. To simulate the effect of varying R.A. in the actual data, MC events are randomly assigned R.A. values under the assumption of a flat local sidereal time. Several corrections are applied on an event-by-event basis to account for neutrino oscillations and to reflect best-fit values of systematic error parameters from the analysis in Abe et al. (2018). The event rates in the search regions are then calculated for each SK phase and combined with appropriate livetime normalization factors.

Figure 2 shows the observed data events in each fixed search cone superimposed on the MC taken over various declinations. Note that the FC and PC background events distribute across all declinations and show a slight increase at higher declinations, especially PC events, due to the decrease of UPMU neutrinos lost to neutrino oscillations. The double-bump structure around −50° and 50° is due to the increased atmospheric neutrino flux coming from the near-horizon direction, where the effective atmospheric depth is deeper than other directions. Since UPMU events are required to come from below the horizon, their maximum decl. is about 54°. The observed data agree with the expected background within 0.7σ for FC, 1.1σ for PC, and 1.2σ for UPMU events, considering statistical uncertainties alone.

Zoom In Zoom Out Reset image size

To further check the consistency of the observed event rates in the search cones, we further investigate by studying similarly sized angular regions taken at the same decl. as TXS 0506+056 but with different R.A. values. The average and variance of the number of observed events in these "off-source" regions are compared with those in the "on-source" region around the blazar. They are consistent within 0.5σ for FC, 1.6σ for PC, and 1.5σ for UPMU based on counting statistics only. The event rate in the on-source search cone is 3.0 ± 0.7 for FC, 4.9 ± 0.9 for PC, and 3.4 ± 0.8 for UPMU events per 1000 livetime days. Averaging the off-source rates yields 2.7 ± 0.6 for FC, 3.9 ± 0.6 for PC, and 2.5 ± 0.6 events per 1000 livetime days. Therefore, the on-source and off-source rates are consistent, indicating no excess of neutrino events in the direction of the blazar.

We additionally searched for evidence of a local increase in the neutrino event rate in the period 1996 April to 2018 February to test for a possible correlation with gamma-ray flaring of the blazar. Since the atmospheric neutrino rate is known to be stable at each SK phase, the number of observed neutrino events is expected to increase linearly with increasing livetime if there is no neutrino emission from the blazar. On the other hand, the event rate would deviate from linearity if there were additional neutrinos from the gamma-ray flare. In order to test for the presence of such a variation, we evaluated the probability (p-value) that the observed rate is consistent using a Kolmogorov–Smirnov test (KS-test). Figure 3 compares the cumulative observed events with the expected events as a function of livetime day. To estimate the degree of deviation, a set of 10,000 pseudo-experiments was generated assuming that the expected background from the MC in each SK phase was distributed according to a Poisson function. In each pseudo-experiment, the observed time of each event was randomly assigned assuming a constant rate in each SK phase. The maximum distance between each pseudo-experiment and the expectation is compared to that from the observed SK data to calculate a p-value. This represents the percentage of pseudo-experiments that have a maximum distance larger than the data. The estimated p-values are 91.74%, 12.26%, and 48.75% for FC, PC, and UPMU events, respectively, indicating consistency with a constant event rate. Accordingly, we conclude that no significant signal from the direction of blazar TXS 0506+056 exists in the SK data during the observation period considered here.

Zoom In Zoom Out Reset image size

Since no significant indication of a signal from the blazar was found in any of the tests above, we estimate flux limits based on the expected background throughout the entire observation period. In the following we derive limits on the neutrino fluence from this blazar. Based on the assumption of no signal, the upper limit on the neutrino fluence is estimated following Swanson et al. (2006) and Thrane et al. (2009):

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{{\rm{FC}}}^{{\nu }_{x}+{\bar{\nu }}_{x}}\\ =\,\displaystyle \frac{{N}_{90}^{{\rm{FC}}}}{{N}_{T}\int {{dE}}_{\nu }\left({\sigma }^{{\nu }_{x}}({E}_{\nu }){\varepsilon }^{{\nu }_{x}}({E}_{\nu })+{\sigma }^{{\bar{\nu }}_{x}}({E}_{\nu }){\varepsilon }^{{\bar{\nu }}_{x}}({E}_{\nu })\right)\lambda ({E}_{\nu }^{-2})}\\ \,\,(x=e,\mu ),\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{PC}}^{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}\\ =\,\displaystyle \frac{{N}_{90}^{\mathrm{PC}}}{{N}_{T}\int {{dE}}_{\nu }\left({\sigma }^{{\nu }_{\mu }}({E}_{\nu }){\varepsilon }^{{\nu }_{\mu }}({E}_{\nu })+{\sigma }^{{\bar{\nu }}_{\mu }}({E}_{\nu }){\varepsilon }^{{\bar{\nu }}_{\mu }}({E}_{\nu })\right)\lambda ({E}_{\nu }^{-2})},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{UPMU}}^{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}\\ =\,\displaystyle \frac{{N}_{90}^{\mathrm{UPMU}}}{{A}_{\mathrm{eff}}(z)\int {{dE}}_{\nu }\left({P}^{{\nu }_{\mu }}({E}_{\nu }){S}^{{\nu }_{\mu }}(z,{E}_{\nu })+{P}^{{\bar{\nu }}_{\mu }}({E}_{\nu }){S}^{{\bar{\nu }}_{\mu }}(z,{E}_{\nu })\right)\lambda ({E}_{\nu }^{-2})}.\end{array}\end{eqnarray} \tag{ 3 }$$

Here N₉₀ is the upper limit on the number of the events above the background at the 90% confidence level (C.L.), N_T is the number of nucleons in the 22.5 kiloton fiducial volume of the detector, σ(E_ν) is the total neutrino interaction cross-section from the NEUT model, and ε(E_ν) is the neutrino detection efficiency. The parameter A_eff is the effective area of SK to UPMU interactions, and P(E_ν) is the probability for a neutrino to produce a muon in the rock surrounding that which reaches SK. Calculation of the latter is done using the charged-current cross-section for muon neutrino–nucleon interaction in rock coupled with the expected range of the resulting muon produced assuming the initial vertex is some distance from the detector. The attenuation of neutrinos due to their interactions in the Earth is given by S(z, E_ν) and taken from Gandhi et al. (1996) using the "Preliminary Earth Model" (Dziewonski & Anderson 1981). Finally, $\lambda ({E}_{\nu }^{-2})$ is the number density distribution from the blazar's direction, which we assume to follow a power law with spectral index of −2 as in the IceCube analysis (Aartsen et al. 2018a).

Assuming Poisson statistics for the expected backgrounds described above, N₉₀ for the FC, PC, and UPMU samples has been estimated to be 10.2, 14.6, and 12.7, respectively. For the FC sample, we conservatively make no distinction between electron- and muon-type neutrinos when calculating N₉₀.⁵⁴ The energy ranges used in the integrals in Equation (3) are 5.1–10 GeV (FC), 1.8–100 GeV (PC), and 1.6 GeV–10 TeV (UPMU). These ranges represent the MC neutrino energies populating each sample. Upper limits are calculated for both electron and muon neutrinos using FC events since this sample is sensitive to both. Only that of muon-neutrino fluence limits are estimated for the other samples. The results are shown in Table 1. In these calculations, the fluence limits are calculated using the average zenith angle to the source taken over the detector observation period, since the effective area, A_eff(z), and the shadowing effect, S(z, E_ν), depend on the zenith angle and fluctuate with the motion of the Earth.

Table 1.  Summary of the Data Samples and Fluence Limits

	FC	PC	UPMU
Energy Range (GeV)	5.1–10	1.8–100	1.6–104
Search Cone (°)	10	10	5
Observed Events	18	29	20
Expected Background	15.2	22.9	14.5
N90	10.2	14.6	12.7
Fluence Limit (${\nu }_{\mu }+{\bar{\nu }}_{\mu }$) (cm−2)	6.9 × 104	1.1 × 105	1.3 × 102
Fluence Limit (${\nu }_{e}+{\bar{\nu }}_{e}$) (cm−2)	1.9 × 104	⋯	⋯

Note. Upper limits at 90% C.L. on the neutrino fluence have been calculated assuming an E⁻² energy spectrum from blazar TXS 0506+056.

Download table as:  ASCIITypeset image

Figure 4 shows the energy-dependent upper limits for electron-neutrino (${\nu }_{e}+{\bar{\nu }}_{e}$) and muon-neutrino (${\nu }_{\mu }+{\bar{\nu }}_{\mu }$) fluxes by SK observations compared by the IceCube group (Aartsen et al. 2018a). The IceCube collaboration considered two neutrino emission periods to calculate the flux limit. In the first scenario, neutrinos were assumed to be emitted only during the about 6 month period corresponding to the duration of the γ-ray flare. Alternatively, neutrinos emitted over the whole observation of IceCube (7.5 yr) were considered. These two benchmark cases and the result of our analysis are shown in Figure 4.

Zoom In Zoom Out Reset image size

Above 10 GeV the upper limits are obtained from UPMU data by using the same formula described above with the corresponding energy range of the integration. It should be noted that the neutrino energy cannot be reconstructed for UPMU events since they are produced by neutrinos interacting in the rock surrounding the detector. Therefore, the same value of N₉₀ is used to calculate the upper limit in each energy bin. Note that the FC sample is populated almost entirely by events with energy less than 10 GeV, so our limit for electron neutrinos spans a single bin. For muon neutrinos in the 1–10 GeV bin and in the 10–100 GeV bin, the limit is calculated using the total observation and expectation, summing over the FC, PC, and UPMU contributions. The flux limit by SK covers ∼1 GeV–10 TeV, and becomes more stringent as the energy increase because of larger neutrino cross-section and smaller atmospheric neutrino backgrounds.

We note that the IceCube group observed evidence of a neutrino event excess between 2014 and 2015 from the direction of the blazar whose best-fit energy spectrum was E^−2.2±0.2 and whose flux was ${2.5}_{-1.0}^{+1.1}\times {10}^{-9}\,[{\rm{e}}{\rm{r}}{\rm{g}}\,{{\rm{c}}{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}]$ at 100 TeV. This corresponds to a flux of 2.0 × 10⁻⁸ [erg cm⁻² s⁻¹] at 3.2 GeV and 5.1 × 10⁻⁹ [erg cm⁻² s⁻¹] at 3.2 TeV, respectively. They can be compared limits from the present analysis of 4.5 × 10⁻⁷ [erg cm⁻² s⁻¹] and 9.3 × 10⁻¹⁰ [erg cm⁻² s⁻¹] in the respective energy regions. If the IceCube spectral index is used in our calculations, the flux limits increase by 15% at most.

We performed a search for neutrino detections coincident in the direction of blazar TXS 0506+056 using the observation data from 1996 April to 2018 February by Super-Kamiokande in GeV to several TeV regions. By comparing to the expected backgrounds, no significant excess was observed at greater than the 2σ level in the blazar direction of 5° (10°) for UPMU (FC, PC) samples. No significant temporal increase of neutrino flux was found in the blazar direction by examining the change of the event rate using the Kolmogorov–Smirnov test. Based on no signal assumption, upper limits of neutrino fluence and the energy-dependent neutrino flux are given for both electron neutrinos and muon neutrinos. Upper limits are placed on the electron-neutrino flux of 6.0 × 10⁻⁷ [erg cm⁻² s⁻¹] below 10 GeV and on the muon-neutrino flux of 4.5 × 10⁻⁷ [erg cm⁻² s⁻¹] to 9.3 × 10⁻¹⁰ [erg cm⁻² s⁻¹] in the range 1 GeV to 10 TeV assuming an E⁻² energy spectrum.

We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super-Kamiokande experiment has been built and operated from funding by the Japanese Ministry of Education, Culture, Sports, Science and Technology, the U.S. Department of Energy, and the U.S. National Science Foundation. Some of us have been supported by funds from the National Research Foundation of Korea NRF-2009-0083526 (KNRC) funded by the Ministry of Science, ICT, and Future Planning and the the Ministry of Education (2018R1D1A3B07050696, 2018R1D1A1B07049158), the Japan Society for the Promotion of Science, the National Natural Science Foundation of China under grants No. 11235006, the National Science and Engineering Research Council (NSERC) of Canada, the Scinet and Westgrid consortia of Compute Canada, the National Science Centre, Poland (2015/18/E/ST2/00758), the Science and Technology Facilities Council (STFC) and GridPPP, UK, and the European Unions H2020-MSCA-RISE-2018 JENNIFER2 grant agreement No. 822070.