MEASUREMENT OF THE ANISOTROPY OF COSMIC-RAY ARRIVAL DIRECTIONS WITH ICECUBE

Long-term observations of cosmic-ray muons by underground experiments have demonstrated the presence of an anisotropy in the cosmic-ray intensity up to a few hundred GeV (Nagashima et al. 1998). Recent underground and surface array measurements of cosmic rays by the Tibet array (Amenomori et al. 2006), Super-Kamiokande (Guillian et al. 2007), and Milagro (Abdo et al. 2009) indicate that the anisotropy persists into the TeV range.

All of the TeV measurements were performed in the Northern Hemisphere; so far, no such measurement has been performed covering the entire Southern Hemisphere at median energies in the multi-TeV region. With the deployment of the IceCube Neutrino Observatory at the South Pole, we have for the first time measured the anisotropy at TeV energies in the southern sky. IceCube is primarily a neutrino detector, but it is sensitive to the muons produced in downward-going cosmic-ray air showers. The observatory provides high-statistics measurements of cosmic rays with median energy of 20 TeV.

When completed in 2011, IceCube will comprise 5160 optical modules buried 1450 and 2450 m below the surface of the polar ice sheet. The modules are physically connected to the surface by electronic umbilical lines, or "strings," with 86 strings in total (Abbasi et al. 2009). In this Letter, we use cosmic-ray data recorded by the detector in its 22 strings configuration (IC22) between 2007 June and 2008 March to produce the cosmic-ray skymap of the Southern sky in the TeV range.

During the IC22 physics run, cosmic-ray events were observed at an average trigger rate of about 550 Hz. The arrival direction is determined by a likelihood-based reconstruction which is seeded with a fast online estimate of the arrival direction (Ahrens et al. 2004). The likelihood-based reconstruction is applied if twelve or more optical sensors on at least three strings were triggered by the event. A total of 5.2 × 10⁹ events satisfied the above conditions at an average rate of ∼240 Hz. Further selection criteria were applied to the data to ensure good quality and stable runs. The final data set contains 4.3 × 10⁹ events with a total livetime of 226 days, a median angular resolution of 3°, and a median energy per cosmic ray of 20 TeV. The energy scale was determined with a standard cosmic-ray simulation program, CORSIKA⁴³, using the SIBYLL hadronic interaction model (Version 2.1; Ahn et al. 2009) and the Poly-Gonato model for the composition and spectrum of the primary cosmic rays (Hörandel 2003).

To evaluate physical anisotropies in the cosmic-ray data set, it is necessary to eliminate spurious effects which can mimic an anisotropy. These include local effects such as diurnal and seasonal variations of atmospheric conditions, asymmetries in the detector geometry, and non-uniform detector exposure to different regions of the sky. Fortunately, the location of IceCube at the South Pole is ideal to compensate for many effects that can impact cosmic-ray detectors in the middle latitudes. At the South Pole, the southern celestial sky is fully visible at any given time, providing complete and uniform coverage. While the seasonal variation in the cosmic-ray event rate is on the order of ±10% (Tilav et al. 2010), these variations are sufficiently slow to have no effect on the anisotropy. Rapid atmospheric changes, which can affect the rate, are rare and can be identified from the data.

The remaining effects, which must be accounted for in this analysis, are an asymmetry in the IceCube detector response, and a non-uniformity in the time coverage of the data. The asymmetric response is due to the geometrical configuration of IceCube during the IC22 physics run; events arriving along the long axis of the detector were preferentially selected by the online filter and reconstruction due to the larger number of strings and modules triggered. In principle, the rotation of the Earth should average out the local asymmetry in the arrival directions each day, but gaps in the detector uptime and uneven run selection due to quality selection introduce non-uniformities into the time coverage of the data. These non-uniformities preclude the complete averaging and translate into an artificial arrival direction asymmetry in equatorial coordinates.

To correct for this detector-related asymmetry, each event from a given local azimuth bin i was weighted with the ratio $\bar{n}/n_i$, where $\bar{n}$ is the average number of events over the full range of local azimuths, and n_i is the number of events in local azimuth bin i. Since the local azimuth distribution varies with zenith angle, the events were grouped into four zenith bands with approximately equal numbers of events per band. The weighting is applied within each band to remove the detector asymmetry.

To investigate the arrival direction distribution of the cosmic rays, we studied the relative intensity of the cosmic-ray-induced muon flux. The arrival direction distribution is dominated by the zenith angle dependence of the muon flux. The zenith angle dependence is a result of varying overburden for the muons through the Antarctic ice, the maximum zenith angle of the events used in this analysis is ∼55°. Due to this zenith dependence, the flux was normalized to unity within declination belts of width 3°, which corresponds to the angular resolution of the data. Note that since the declination belts in the equatorial map are treated independently, the map provides only information on the relative modulation of the arrival direction of cosmic rays along the right ascension.

Figure 1 shows the relative intensity of the event rate in equatorial coordinates. The color scale quantifies the number of reconstructed events with respect to the average number of events in each declination belt. Figure 2 shows the same data in Galactic coordinates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows an anisotropy that appears to be a continuation of a similar modulation of the cosmic-ray flux observed in the Northern Hemisphere (Amenomori et al. 2006; Guillian et al. 2007; Abdo et al. 2009). To quantify the scale of the anisotropy, we fitted the declination belts of the data to a first- and second-order harmonic function of the form

$$\begin{equation} \sum _{i=1}^{n=2}A_{i}\cos (i(\alpha -\phi _{i})) + B, \end{equation} \tag{ 1 }$$

where (A_i, ϕ_i) are the amplitude and phase of the anisotropy, α is the right ascension, and B is a constant. The results of the first- and second-harmonic fit of Equation (1) for the declination belts are given in Table 1. Moreover, the two-dimensional relative intensity skymap is projected to one-dimensional projection in right ascension as shown in Figure 3. The profile in right ascension is obtained by calculating the average value of the relative intensity for the bins that lie within each of the right ascension bands. Each right ascension band is represented in a point in Figure 3. The error bars are derived by propagating the statistical errors from the bins in the right ascension band, and the gray-shaded area indicates the estimated spread from the fit values of the stability tests. The solid line indicates the fit of Equation (1) to the data. The first- and second-harmonic fit parameters to the one-dimensional projection in Figure 3 are A₁ = (6.4 ± 0.2 stat. ± 0.8 syst.) × 10⁻⁴, ϕ₁ = 664 ± 26 stat. ± 38 syst., A₂ = (2.1 ± 0.3 stat. ± 0.5 syst.) × 10⁻⁴, ϕ₂ = −656 ± 40 stat. ± 75 syst., and B = (1 ±  2.1) × 10⁻⁵ stat. ±  9.9 × 10⁻⁶ syst. with χ²/dof = 22/19.

Zoom In Zoom Out Reset image size

To estimate the stability of the result and the corresponding systematic uncertainties, we performed two types of tests: we checked whether the observation is stable against the choice of the particular event sample selection, and whether the modulation has spurious influences from other physical effects. The stability of the result was tested through a series of dedicated checks. The first stability test was done by dividing the data set in half by sub-run number, where each sub-run contains approximately 20 minutes of observations. To avoid any systematic biases, the division was tried in several ways: first, by separating data in even- and odd-numbered sub-runs; and second, by random selection of half of the sub-runs. The corresponding relative intensity distributions in right ascension for both tests were determined, and it was found that the variations induced by the data set selection are within the statistical fluctuations. In addition, to check for daily variational effects, the data were divided in two sets: the first containing sub-runs with event rates above the median value for the corresponding day, and the second containing sub-runs with event rate below the median value. The corresponding relative intensity modulations were fit and found to be smaller than the statistical fluctuations. This means that a variation in the absolute event rate does not affect the modulation on arrival direction in right ascension.

Similarly, more stability tests were applied to check for effects due to seasonal variations and time gaps. To check for the seasonal effect, data were divided into one set containing the winter months (June–October) and the other set containing the summer months (November–March). The relative intensity variations in right ascension were fit and found to be consistent with the statistical fluctuations. To verify that the non-uniform time coverage due to missing sub-runs and other gaps in the data is correctly handled by the azimuthal re-weighting procedure (see Section 2), the relative intensity distribution from the full data set was compared with the one determined using only the days with minimal time gaps. The differences were found to be consistent with statistical fluctuations.

In each of the above stability tests, an independent fit was made to the relative intensity distribution as a function of right ascension using Equation (1). The envelope from all the stability tests' fit curves was constructed, and it is shown as the gray-shaded area in Figure 3.

To verify whether the analysis procedure could induce a modulation in right ascension, the experimental event arrival directions were randomized to generate an isotropic distribution, and the same analysis was performed on this sample. The result was found to be consistent with isotropy.

To check whether the observed anisotropy has some sidereal spurious effect derived from the interference between possible yearly modulated daily variations, the same analysis was performed using the anti-sidereal time frame (a non-physical time defined by switching the sign of the transformation from universal to sidereal time; Farley & Storey 1954). The real feature in the sidereal time is expected to be scrambled in the anti-sidereal time. Figure 4 shows the one-dimensional projection in right ascension for the sidereal time in black and for the anti-sidereal time in red. The amplitude of the first-harmonic fit to the one-dimensional projection in the anti-sidereal time was found to be 0.8 × 10⁻⁴. This value is larger than the spread found in the first-harmonic amplitude from the stability tests, therefore we use it as the systematic uncertainty in the first-harmonic amplitude. The uncertainty in the first-harmonic phase implied by the study in the anti-sidereal time frame is within the systematic error determined from the stability tests. The systematic uncertainties for the rest of the parameters of the fit (quoted in Section 3) are derived from the stability tests.

Zoom In Zoom Out Reset image size

Using a high-statistics sample of downgoing cosmic rays with a median energy of 20 TeV and a median angular resolution of 3°, we presented the first map of the relative intensity of the flux of TeV cosmic rays in the southern sky. The arrival direction distribution of the cosmic rays is found to be anisotropic with a first-harmonic amplitude and phase of A₁ = (6.4 ± 0.2 stat. ± 0.8 syst.) × 10⁻⁴ and ϕ₁ = 664 ± 26 stat. ±  38 syst. The observation appears to be a continuation of a previously measured cosmic-ray anisotropy reported in the Northern Hemisphere (Amenomori et al. 2006; Guillian et al. 2007; Abdo et al. 2009).

The origin of the anisotropy remains unclear. The Compton–Getting effect (Compton & Getting 1935) suggests that the relative motion of the solar system around the Galactic center in the cosmic-ray plasma should give rise to an excess in the direction of motion of the solar system and a deficit in the opposite direction. In this model, an excess flux should appear with a maximum in right ascension between 290° and 340° and a minimum in right ascension between 110° and 160° (Amenomori et al. 2006). As shown in Figure 4, the excess cannot be described in terms of the direction of motion of the solar system. Therefore, we conclude that the Compton–Getting effect could be (at most) one of several contributions to the cosmic-ray anisotropy. This effect will be addressed in more detail in a future study including the energy dependence of the anisotropy.

It is tempting to try to interpret the cosmic-ray excess as an artifact of the heliospheric magnetic field. However, the maximum gyroradius of a 10 TeV cosmic-ray proton in a 1 μG magnetic field is about 0.01 pc, i.e., much larger than the size of the heliosphere. As a consequence, the observed anisotropy is more likely to be connected to features of the local interstellar magnetic field at distances <1 pc. We are also investigating the possibility that the cosmic-ray excess is associated with structures in the Galactic magnetic field at larger distance scales, or with diffusive particle flows from a nearby Galactic source such as Vela.

The still growing IceCube observatory will be completed in 2011 with a total of 86 strings and a volume of 1 km³. The estimated rate of cosmic-ray-induced muons will be greater than 30 billion events per year. Such high statistical power, together with an estimated energy resolution of about 0.3 in log(E), will allow us to determine, in one year, the variation of cosmic-ray anisotropy in several energy ranges up to a few hundred TeV. The energy dependence study will provide fundamental hint at the nature of the source or sources of the cosmic rays, as well as their propagation through the Galactic magnetic field.

We acknowledge the support from the following agencies: U.S. National Science Foundation—Office of Polar Program, U.S. National Science Foundation—Physics Division, University of Wisconsin Alumni Research Foundation, U.S. Department of Energy, and National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; 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), 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); Marsden Fund, New Zealand; Japan Society for Promotion of Science (JSPS); the Swiss National Science Foundation (SNSF), Switzerland; A. Groß acknowledges support by the EU Marie Curie OIF Program; J. P. Rodrigues acknowledges support by the Capes Foundation, Ministry of Education of Brazil.