INDICATIONS OF INTERMEDIATE-SCALE ANISOTROPY OF COSMIC RAYS WITH ENERGY GREATER THAN 57 EeV IN THE NORTHERN SKY MEASURED WITH THE SURFACE DETECTOR OF THE TELESCOPE ARRAY EXPERIMENT

The origin of ultrahigh-energy cosmic rays (UHECRs), particles with energies greater than 10¹⁸ eV, is one of the mysteries of astroparticle physics. Greisen, Zatsepin, and Kuz'min (GZK) predicted that UHECR protons with energies greater than ∼60 EeV (6 × 10¹⁹ eV) would be severely attenuated primarily due to pion photoproduction interactions with the cosmic microwave background radiation (Greisen 1966; Zatsepin & Kuz'min 1966). This GZK suppression becomes strong if these very high energy cosmic rays are produced at and traveling moderate extragalactic distances. The High Resolution Fly's Eye (HiRes) collaboration was first to observe a suppression of cosmic rays above ∼60 EeV Abbasi et al. (2008), which is consistent with expectation from the GZK cutoff. This suppression was independently confirmed by both the Pierre Auger Observatory (PAO; Abraham et al. 2008) in the south and Telescope Array (TA) experiment (Abu-Zayyad et al. 2013a) in the north, which are the largest aperture cosmic-ray detectors currently in operation.

The distribution of UHECR sources should be limited within the local universe with distances smaller than 100 Mpc for proton/iron and 20 Mpc for helium/carbon/nitrogen/oxygen (distances within which ∼50% of cosmic rays are estimated to survive; Kotera & Olinto 2011). To accelerate particles up to the ultrahigh-energy region, particles must be confined to the accelerator site for more than a million years by a magnetic field and/or a large-scale confinement volume (Hillas 1984; Ptitsyna & Troitsky 2010). This would thus limit the number of possible accelerators in the universe to astrophysical candidates such as galaxy clusters, supermassive black holes in active galactic nuclei (AGNs), jets and lobes of active galaxies, starburst galaxies, gamma-ray bursts, and magnetars. Galactic objects are not likely to be the sources since past observations indicate that the UHECRs do not concentrate in the galactic plane and have a relatively isotropic distribution. In addition, our Galaxy cannot confine UHECRs above 10¹⁹ eV within its volume by the GMF. Extragalactic astrophysical objects form well-known large-scale structures (LSSs), most of which are spread along the "supergalactic plane" in the local universe. Nearby AGNs are clustered and concentrated around LSS with a typical clustering length of 5–15 Mpc, as observed by Swift BAT (Cappelluti et al. 2010). Concentrations of nearby AGNs coincide spatially with the LSS of matter in the local universe, including galaxy clusters such as Centaurus and Virgo. The typical amplitude of such AGN concentrations is estimated to be a few hundred percent of the averaged density within a 20° radius circle, which is of an angular scale comparable to the clustering length of the AGNs within 85 Mpc (Ajello et al. 2012).

The main difficulty in identifying the origin of UHECRs is the loss of directional information due to magnetic field induced bending. In order to investigate the UHECR propagation from the extragalactic sources, a number of numerical simulations have been developed (Yoshiguchi et al. 2003; Sigl et al. 2004; Takami et al. 2006; Kashti & Waxman 2008; Koers & Tinyakov 2009; Takami & Sato 2010; Kalli et al. 2011; Takami et al. 2012). In the simulations, the UHECR trajectory between the assumed UHECR source and the Earth is traced through intergalactic and galactic magnetic fields (IGMF and GMF). The results depend strongly on the assumed distribution and density of the UHECR sources and the intervening magnetic fields. The deflection angle of a 60 EeV proton from a source at a distance of 50 Mpc is estimated to be a few degrees assuming models with an IGMF strength of 1 nG. Meanwhile, the estimated deflection by the GMF ranges from a few to about 10°. This, however, depends on the direction in the sky. If the highest-energy cosmic rays come from the local universe such as nearby galaxies, and if they are protons, the maximum amplitude of the cosmic-ray anisotropy above ∼60 EeV is expected to be a few hundred percent of the average cosmic-ray flux. In this case, the amplitude of the cosmic-ray anisotropy might be detectable by the UHECR detectors of the TA and PAO.

In the highest-energy region, E > 57 EeV, the PAO found correlations of the cosmic-ray directions within a 31 radius circle centered at nearby AGNs (within 75 Mpc) in the southern sky (Abraham et al. 2007). Updated measurements from the PAO indicate a weakened correlation with nearby AGNs (Abreu et al. 2010; Macolino 2012); the correlating fraction (the number of correlated events divided by all events) decreased from the early estimate of ($69^{+11}_{-13}$)% to (33 ± 5)%, compared with 21% expected for an isotropic distribution of cosmic rays. The chance probability of the original (69%) correlation is 6 × 10⁻³ assuming an isotropic sky. The TA has also searched for UHECR anisotropies such as autocorrelations, correlations with AGNs, and correlations with the LSS of the universe using the first 40 months of scintillator surface detector (SD) data (Abu-Zayyad et al. 2012b, 2013b). Using 5 yr of SD data, we updated results of the cosmic-ray anisotropy with E > 57 EeV, which shows deviations from isotropy at the significance of 2σ–3σ (Fukushima et al. 2013). In this Letter, we report on indications of intermediate-scale anisotropy of cosmic rays with E > 57 EeV in the northern hemisphere sky using the 5 yr TA SD data set.

The TA is the largest cosmic-ray detector in the northern hemisphere. It consists of a scintillator SD array (Abu-Zayyad et al. 2012a) and three fluorescence detector (FD) stations (Tokuno et al. 2012). The observatory has been in full operation in Millard Country, Utah, USA (3930N, 11291W; about 1400 m above sea level) since 2008. The TA SD array consists of 507 plastic scintillation detectors each 3 m² in area and located on a 1.2 km square grid. The array has an area of ∼700 km². The TA SD array observes cosmic-ray-induced extensive air showers with E > ∼1 EeV, regardless of weather conditions with a duty cycle near 100% and a wide field of view (FoV). These capabilities ensure a very stable and large geometrical exposure over the northern sky survey in comparison with FD observations that have a duty cycle of ∼10%.

In this analysis, we used SD data recorded between 2008 May 11 and 2013 May 4. The dataset contains approximately 1 million triggered events. For the reconstructed events, the energies determined by the SD array were renormalized by 1/1.27 to match the SD energy scale to that of the FD, which was determined calorimetrically (Abu-Zayyad et al. 2013a). Of these events, 72 met the following conditions: (1) each event included at least four SD counters, (2) the zenith angle of the event arrival direction was less than 55°, and (3) the reconstructed energy was greater than 57 EeV, which corresponds to the energy threshold determined from the AGN correlation analysis results obtained by the PAO (Abraham et al. 2007), and is adopted here to avoid introducing a free parameter in the scanning phase space.

The event selection criteria above are somewhat looser than those of our previous analyses of cosmic-ray anisotropy (Fukushima et al. 2013) to increase the observed cosmic-ray statistics. In our previous analyses, the largest signal counter is surrounded by four working counters that are its nearest neighbors to maintain the quality of the energy resolution and angular resolution. Only 52 events survived those tighter cuts. When the edge cut is abolished from the analysis (presented here) to keep more cosmic-ray events, 20 events with E > 57 EeV are recovered compared with the tighter cut analysis. A full Monte Carlo (MC) simulation, which includes detailed detector responses (Abu-Zayyad et al. 2013a), predicted a 13.2 event increase in the number of events. The chance probability of the data increment being 20 as compared to the MC prediction of 13.2 is estimated to be 5%, which is within the range of statistical fluctuations. The angular resolution of array boundary events deteriorates to 17, compared to 10 for the well contained events. The energy resolution of array boundary events also deteriorates to ∼20%, where that of the inner array events is ∼15%. These resolutions are still good enough to search for intermediate-scale cosmic-ray anisotropy. One final check is that when we calculate the cosmic-ray spectrum using the loose cuts analysis, the result is consistent with our published spectrum.

Figure 1(a) shows a sky map in equatorial coordinates of the 72 cosmic-ray events with energy E > 57 EeV observed by the TA SD array. A cluster of events appears in this map centered near right ascension ∼150°, and declination ∼40°, with a diameter of ∼30°–40°. In order to determine the characteristics of the cluster, and estimate the significance of this effect, we choose to apply elements of an analysis that was developed by the AGASA collaboration to search for large-size anisotropy (Hayashida et al. 1999a, 1999b), namely to use oversampling with a 20° radius. Being mindful that scanning the parameter space of the analysis causes a large increase in chance corrections, we have not varied this radius. The TA and HiRes collaborations used this method previously (Kawata et al. 2013; Ivanov et al. 2007) to test the AGASA intermediate-scale anisotropy results with their data in the 10¹⁸ eV range. The present letter reports on an extension of this method with application to the E > 57 EeV energy region.

Zoom In Zoom Out Reset image size

In our analysis, at each point in the sky map, cosmic-ray events are summed over a 20° radius circle as shown in Figure 1(b). The centers of tested directions are on a 01 × 01 grid from 0° to 360° in right ascension (R.A.) and −10°–90° in declination (decl.). We found that the maximum of N_on, the number of observed events in a circle of 20° radius is 19 within the TA FoV. To estimate the number of background events under the signal in N_on, we generated 100,000 events assuming an isotropic flux. We used a geometrical exposure g(θ) = sin θcos θ as a function of zenith angle (θ) because the detection efficiency above 57 EeV is ∼100%. The zenith angle distribution deduced from the geometrical exposure is consistent with that found in a full MC simulation. The MC generated events are summed over each 20° radius circle in the same manner as the data analysis, and the number of events in each circle is defined as N_off. Figure 1(c) shows the number of background events N_bg = ηN_off, where η = 72/100,000 is the normalization factor.

We calculated the statistical significance of the excess of events compared to the background events at each grid point of sky using the following equation (Li & Ma 1983):

$$\begin{equation} S_{\rm LM} \!=\! \sqrt{2}\left[N_{\rm on} {\rm ln} \! \left(\frac{(1+\eta)N_{\rm on} }{\eta (N_{\rm on} + N_{\rm off})}\right) + N_{\rm off} {\rm ln} \! \left(\frac{(1+\eta) N_{\rm off} }{N_{\rm on} + N_{\rm off}} \right)\right]^{1/2}\!. \end{equation} \tag{ 1 }$$

Figure 1(d) shows a significance map (in equatorial coordinates) of the events above 57 EeV as observed by the TA SD array. The maximum excess in our FoV appears as a "hotspot" centered at R.A.(α) =1467, decl.(δ) =432 with a statistical significance of S_MAX = 5.1σ (N_on = 19, N_bg = 4.49).

The significance of the hotspot, quoted above at 5.1σ, does not take random clustering into account, so one must make a correction. We did not carry out a blind analysis, but have been watching the hotspot grow over several years as we collected further data and added events to the sky plot. It is difficult to estimate the penalty due to our having seen the cluster of events. For example, in applying the oversampling technique used by the AGASA experiment, we knew the oversampling radius roughly matched the size of the hotspot cluster.

However, by making a simple MC calculation, one can estimate the probability of such a hotspot appearing by chance anywhere in an isotropic sky. One generates many isotropic MC event sets, each with the statistics of the experimental data, then performs a calculation of the Li–Ma significance exactly as was done on the data; i.e., using oversampling with a radius of 20°. One can go further and approximate the effect of the eye's estimate of the radius of the cluster of events by repeating the calculation at other oversampling radii. We did this, choosing five oversampling radii, 15°, 20°, 25°, 30°, and 35°. We chose a 5° scan since by eye one cannot make an estimate more accurately than about ±5°.

We generated 1 million MC data sets, each having 72 spatially random events within our FOV (i.e., we reproduced the statistics of the experimental data), assuming a uniform distribution over the TA SD exposure. The maximum of the significances, S_MAX, was calculated for each MC data set in the same way as in the data, with the exception that the five oversampling radii were used, and the largest S_MAX was chosen. The distribution of the largest S_MAX of the 1 million data sets is shown in Figure 2. We found that there were 365 instances of S_MAX > 5.1σ. This yields a chance probability of the observed hotspot in an isotropic cosmic-ray sky of 3.7 × 10⁻⁴, equivalent to a one-sided probability of 3.4σ.

Zoom In Zoom Out Reset image size

To estimate the size of the hotspot, we present (see Figure 3) the normalized number of events as a function of the opening angle, ψ, relative to the center of the hotspot in the data. Although with current statistics we cannot determine the shape of the hotspot, to estimate its overall size we fit the hotspot excess using the binned maximum likelihood method, assuming a Gaussian signal plus a background estimated by the MC simulation. We used the following equation:

$$\begin{equation} f(\psi; A_{\rm s}, \sigma _{\rm s}) = A_{\rm s} {\rm exp}\left(-\frac{\psi ^{2}}{2\sigma _{\rm s}^{2}}\right) + (a_{0} + a_{1}\psi ^{2} + a_{2}\psi ^{4}), \end{equation} \tag{ 2 }$$

where the first term is the Gaussian signal, and A_s and σ_s denote fitting parameters of the signal height and spread, respectively. The second term is the shape of the background fitted by a polynomial function determined from the MC simulations (a₀ = 0.118, a₁ = −1.7 × 10⁻⁵, and a₂ = 8.5 × 10⁻¹⁰). The spread of the hotspot was σ_s = 103 ± 19 (A_s = 0.67 ± 0.29). The uncertainty in the position of the hotspot is estimated to be $\sigma _{\rm s} / \sqrt{N_{\rm on} - N_{\rm bg}} = 2{.\!\!^{\circ} }7$.

Zoom In Zoom Out Reset image size

There are no known specific sources behind the hotspot. The hotspot is located near the supergalactic plane, which contains local galaxy clusters such as the Ursa Major cluster (20 Mpc from Earth), the Coma cluster (90 Mpc), and the Virgo cluster (20 Mpc). The angular distance between the hotspot center and the supergalactic plane in the vicinity of the Ursa Major cluster is ∼19°.

Assuming the hotspot is real, two possible interpretations are that it may be associated with the closest galaxy groups and/or the galaxy filament connecting us with the Virgo cluster (Dolag et al. 2004); or, if cosmic rays are heavy nuclei, they may originate close to the supergalactic plane and be deflected by extragalactic magnetic fields and the galactic halo field (Tinyakov & Tkachev 2002; Takami et al. 2012). To determine the origin of the hotspot, we will need greater UHECR statistics in the northern sky. Better information about the mass composition of the UHECRs, GMF, and IGMF would also be important.

Using cosmic-ray events with energy E > 57 EeV, collected over five years with the TA SD, we have observed a cluster of events, which we call the hotspot, with a statistical significance of 5.1σ (N_on = 19, N_bg = 4.49), centered at R.A. = 1467, ${\rm {\rm d}ec{\rm l}.}=43{.\!\!^{\circ} }2$. We calculated the probability of such a hotspot appearing by chance in an isotropic cosmic-ray sky to be 3.7 × 10⁻⁴ (3.4σ).

This indication of intermediate-scale anisotropy is limited by statistics collected by experiments in the northern hemisphere. It provides a strong impetus for an improved effort to study the origin of UHECRs. The TA×4 project (extension of the TA SD by a factor of 4; Sagawa 2013) is designed to provide the equivalent of 20 TA-years of SD data by 2019, which would yield a ∼7σ observation if the ratio of hotspot to background events remains as is currently seen. TA×4 and other related projects will enable us to make a precise UHECR anisotropy map with high statistics and help solve the mystery of the UHECR origin.

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science through Grants-in-Aids for Scientific Research on Specially Promoted Research (21000002) "Extreme Phenomena in the Universe Explored by Highest Energy Cosmic Rays" and for Scientific Research (19104006), and the Inter-University Research Program of the Institute for Cosmic Ray Research; by the U.S. National Science Foundation awards PHY-0307098, PHY-0601915, PHY-0649681, PHY-0703893, PHY-0758342, PHY-0848320, PHY-1069280, and PHY-1069286; by the National Research Foundation of Korea (2007-0093860, R32-10130, 2012R1A1A2008381, 2013004883); by the Russian Academy of Sciences, RFBR grants 11-02-01528a and 13-02-01311a (INR), IISN project No. 4.4509.10 and Belgian Science Policy under IUAP VII/37 (ULB). The foundations of Dr. Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles and the George S. and Dolores Dore Eccles all helped with generous donations. The State of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management, and the U.S. Air Force. We also wish to thank the people and the officials of Millard County, Utah for their steadfast and warm support. We gratefully acknowledge the contributions from the technical staffs of our home institutions. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged.

In this Appendix we present the list of events with energy E > 57 EeV and zenith angle θ < 55° that have been recorded by the SD of the TA from 2008 May 11 to 2013 May 4. During this period, 72 such events were observed. Table 1 shows the arrival date and time of these events, the zenith angle, energy in units of EeV, and equatorial coordinates α (R.A.) and δ (decl.) in degrees. Noted that the air shower reconstruction used here as described in Kawata et al. (2013) was slightly different from that of previous anisotropy work (Abu-Zayyad et al. 2012b). The opening angles between these arrival directions and previous ones are almost within 1°. This difference hardly affects the results of the large-scale anisotropy.