Discovery of Localized Regions of Excess 10-TeV Cosmic Rays

Phys. Rev. Lett. 101, 221101
Map of significances for the Milagro data set without any cuts to remove the hadronic cosmic-ray background. A 10° bin was used to smooth the data, and the color scale gives... the significance. The solid line marks the Galactic plane, and every 10° in Galactic latitude are shown by the dashed lines. The black dot marks the direction of the heliotail, which is the direction opposite the motion of the solar system with respect to the local interstellar matter. The fractional excess of region A is 6×104, while for region B it is 4×104. The deep deficits bordering the regions of excess appear because the background calculation has been raised by the excess. Show more

The flux of charged cosmic rays at TeV energies is known to be nearly isotropic. This is due to Galactic magnetic fields, which randomize the directions of charged particles. However, numerous experiments across a wide range of energies have found anisotropy on large angular scales, typically with a fractional amplitude of 10-3 (see [1–5], for example). Large-scale anisotropy is also seen in data from the Milagro detector [6]; here we present the results of an analysis sensitive to intermediate angular scales ( 10°).

Milagro [7] is a water Cherenkov air shower detector located in New Mexico, USA at an altitude of 2630 m and at 36° N latitude. It is composed of a central 60m×80m pond surrounded by a sparse 200m×200m array of 175 “outrigger” water tanks. The pond is instrumented with 723 photomultiplier tubes (PMTs) in two layers. The top layer and outrigger tanks are used to determine the direction and energy, while the bottom layer is used to distinguish between gamma-ray induced and hadron induced air showers. The outriggers, with each tank containing a single PMT, improve the angular and energy resolution of the detector for events collected after May, 2003. Milagro has a 2sr field of view, operates with a >90% duty cycle, and has a trigger rate from cosmic rays of 1700Hz, making it well suited to searching for anisotropy in the arrival directions of TeV cosmic rays.

For studies on small to intermediate scales ( 10°), an adaptation of the gamma-ray point source analysis, which has been published previously [7], is used. The primary difference between the previous analysis and the current analysis is that no cosmic-ray background rejection cuts are made. These cuts removed over 90% of the events, so the analysis reported here uses nearly 10 times the number of events of the previous analysis. Like the previous analysis, a signal map is made based on the arrival direction of each event. A background map is also created using the “direct integration” technique [7], in which two-hour intervals are used to calculate the background. Because of this two-hour interval, the analysis is relatively insensitive to features larger than 30° in right ascension (RA); a different analysis of the Milagro data sensitive to larger features has been performed and is presented elsewhere [6].

Signal and background events vs RA for 10°<δ<20°. The plot was made using independent 10° δ by 1° RA bins (i.e., no smoothing). A subset of the... data was used in which there are only full days of data in order to give an approximately uniform exposure in RA. Region A corresponds to the excess at RA70°. This plot shows that the region A excess is inherent in the raw signal data and is not due to an underestimation of the background. Show more

In the gamma-ray point source analysis, the signal and background maps are smoothed with a square bin of size 2.1°/cos(δ) in RA by 2.1° in declination ( δ), which is optimal for Milagro’s angular resolution. However, the bin size may be increased to improve the sensitivity to larger features, with a maximum size of about 10° for δ<60° (for δ>60°, the RA bin width 10°/cos(δ) becomes too large for the 30° background interval). The significance is calculated using the method of Li and Ma [8].

The analysis has been applied to data collected between July 2000 and August 2007. Events were required to have a zenith angle <45° and nFit 20, where nFit is the number of PMTs used in the angle fit. With these cuts, the data set consists of 2.2×1011 events with a median energy of 1TeV and an average angular resolution of <1°. Figure 1 shows the map of significances made with 10° smoothing and no cuts to discriminate gamma rays from charged cosmic rays. The Cygnus region, which has previously been shown to emit TeV gamma rays [9], is clearly visible. The excesses labeled “region A” and “region B” are seen with peak significances of 15.0σ and 12.7σ, respectively. These are pretrial significances because the location and extent of the excesses were determined by examining the data. A map such as shown in Fig. 1 has a few 100 000 independent bins, but given the high statistical significance many maps could be examined and the post-trial’s significance would be reduced by <1σ. The fractional excess relative to the cosmic-ray background is 6×10-4 for region A and 4×10-4 for region B. Note that both excesses are paralleled by regions of deep deficit; this is a known effect of the analysis due to the fact that regions A and B are included in the background calculation of neighboring areas in RA. Therefore, the excess raises the background calculation above its actual value, resulting in an apparent deficit.

Similarity is seen between the map in Fig. 1 and results from the Tibet AS γ collaboration [3,10], but a direct comparison cannot be made because the analysis methods differ. For each band in δ, the Tibet analysis measured the excess (or deficit) relative to the average for that δ band, making it sensitive to the large-scale anisotropy discussed in [3]. Smaller features, such as regions A and B, were superimposed on the large-scale variation, which is several times greater in amplitude. Conversely, in the analysis presented here, the excess or deficit was measured with respect to the local background calculation, which is determined from the data ±30° in RA. This is illustrated in Fig. 2, which shows the data and background calculation versus RA for a 10° band in declination without any smoothing applied to the data. The large-scale variation dominates the figure, but the background calculation makes the analysis sensitive only to features with an extent smaller than 30° in RA. It is noteworthy that the Tibet AS γ collaboration has developed a model for the large-scale structure, and the residual map after subtracting that model from their data shows excesses similar to regions A and B [10].

(a) Differential plot of the fractional excess versus loge(fout) for regions A and B, where fout is the fraction of the outriggers... hit. The spectrum of region A is significantly different than the background (2×106), which is represented by the horizontal line. (b) Profile plot of the simulated energy of protons for the as a function of loge(fout). The ranges are asymmetric and contain the inner 68% of simulated events. Show more

To estimate the extent of region A, an elliptical Gaussian was fit to the excess map of the data in 0.1° bins prior to smoothing. The fit, which accounted for the change in sensitivity with declination, returned a centroid of RA=69.4°±0.7°, δ=13.8°±0.7°, a half width of 2.6°±0.3°, a half length of 7.6°±1.1°, and an angle of 46°±4° with respect to the RA axis. It is important to note that this fit focused on a “hot spot” in the general excess of region A, but there is still excess extending to lower declinations. A fit was not performed to the excess in region B due to its large, irregular shape.

While the excesses in regions A and B are statistically significant, systematic causes must be ruled out. Potential weather-related effects were explored by dividing the data into the four seasons, and both excesses were seen in each season. The data were also divided into yearly data sets to investigate whether changes to the detector could play a role, and again the excesses were found in each data set. The analysis was also run using universal time instead of sidereal time to check for day-night effects which could masquerade as a signal. In addition, the data were analyzed using antisidereal time, which provides a sanity check on the analysis since it will scramble real celestial features. No excess appears in either analysis.

Potential errors in the background calculation were also investigated. Figure 2 shows the number of events versus RA for the signal and background for 10°<δ<20°, using independent 10°δ by 1° RA bins (i.e., no smoothing). The data for this figure were chosen to include only full days in order to achieve an approximately uniform exposure as a function of RA, and the broad deficit seen by the Tibet Air Shower Array is evident (centered around RA=180°). As can be seen, the background estimate as calculated via the direct integration technique [7] agrees well with the data. The excess corresponding to region A is clearly inherent in the raw signal data and is not an artifact created by the background subtraction. A similar result is found for region B.

Results of a χ2 fit to the excesses in regions A and B, assuming a pure-proton spectrum of the form in Eq. (1). The top... panel shows the results for region A and the bottom panel shows the results for region B. The 1σ, 2σ and 3σ allowed regions of the spectral index γ and the cutoff energy Ec are indicated by the shaded regions. Show more

Diagnostic tests have been performed to gain insight into the nature of regions A and B. For the purposes of these tests, region A is defined as the box bounded by 66°<RA<76° and 10°<δ<20°. Region B is defined as the union of two boxes: 117°<RA<131° and 15°<δ<40°, and 131°<RA<141° and 40°<δ<50°. These definitions were chosen by visual inspection of Fig. 1.

To check for flux variation, the analysis was applied to yearly and seasonal data sets. For each region, the yearly excess was consistent with a constant flux. Both regions also had a significant excess during each of the four seasons, with the respective fractional excess in parts per 10 000 in spring, summer, fall, and winter of 6.5±0.9, 4.0±0.9, 6.4±0.9, 7.2±0.9 for region A and 3.5±0.4, 3.3±0.4, 4.0±0.4, 4.7±0.4 for region B. In both cases the fractional excess was lowest in the summer and highest in the winter, and the χ2 probability relative to a constant fractional excess is only about 5% for each region. While this may provide insight into the cause of these excesses, only statistical errors are given. There could be systematic effects such as the slightly higher energy threshold of Milagro in winter when there is snow on top of the pond.

The excesses in regions A and B are inconsistent with pure gamma-ray emission. We can statistically separate gamma-ray events from cosmic-ray events utilizing two parameters. The compactness parameter [7] uses PMT information in the bottom layer of Milagro to identify the penetrating particles characteristic of a hadronic air shower. The distribution of compactness depends on the energy spectrum of the source with higher energy gamma rays producing showers of greater compactness. In order to exclude a gamma-ray hypothesis of any spectrum, we also fit an energy parameter fout, the fraction of outrigger PMTs that detect light. Figure 3(a) shows the fractional excess of loge(fout) for regions A and B. We hypothesize a spectrum for the excess of the form

where γ is the spectral index and Ec is the characteristic energy at which the spectrum cuts off. We attempt this fit for regions A and B assuming separately that the primary particles are purely gamma rays and purely protons. The gamma-ray hypothesis in region A has a χ2 of 124.0 with 16 degrees of freedom. The cosmic-ray hypothesis produces a reasonable fit, with a minimum χ2 of 10.3. For region B, the best gamma-ray hypothesis has a χ2 of 84.8 compared to a best cosmic-ray hypothesis of χ2=19.0, again with 16 degrees of freedom. Thus the proton hypothesis is a reasonable fit for both regions and the gamma-ray hypothesis is inconsistent with probabilities of 9×10-19 for region A and 2×10-11 for region B. The possibility that the regions contain some admixture of protons and gamma rays has not been considered. Figure 4 shows the 1σ, 2σ, and 3σ regions around the best fit for region A and region B for the pure-proton hypothesis. Some care must be taken in interpreting Fig. 4. It does not account for our systematic errors. There is a estimated systematic uncertainty in the spectral index of ±0.2 due to variation in the trigger threshold (caused by such things as changes in atmospheric pressure or ice on the pond). There is also a 30% systematic uncertainty in the energy scale due to the threshold variation, as well as discrepancy between the simulated and measured trigger rates. The fit does not constrain the spectrum well except to suggest that a hard spectrum is favored, particularly for region A. The cutoff energy is constrained, with log10(Ec/GeV)=4.0-0.5(stat)+0.4 for region A and log10(Ec/GeV)=4.0-0.5(stat)+0.3 for region B. Most importantly, the pure gamma-ray hypothesis is strongly disfavored.

We can see that the excesses in regions A and B are harder than the spectrum of the isotropic part of cosmic rays with minimal systematic effects by looking at the data alone. Figure 3(a) shows the fractional excess in regions A and B as a function of fout. Assuming that the excess is due to cosmic rays of the same spectrum, we would expect the fractional excess to be completely flat. The offset from zero tells us that this region is in fact an excess. A χ2 test of whether the curves in Fig. 3 are flat for region A(B) returns a chance probability of 2×10-6 ( 6×10-3), independent of systematic errors. The excess of region A is most significantly detected for loge(fout)-1.5, corresponding to an energy of about 10 TeV for protons, as shown in Fig. 3(b). At around 10 TeV, the spectrum cuts off consistent with the results of the spectral fit.

There is currently no compelling explanation for the excesses in regions A and B. One possibility is that they could be due to neutrons, but this is unlikely because the decay length of 10 TeV neutrons is only about 0.1 parsecs, which is much closer than the nearest star. Another possibility is that these excesses could be caused by a Galactic cosmic-ray accelerator, but this is difficult because the gyroradius of a 10 TeV proton in a 2μG magnetic field, which is the estimated strength of the local Galactic field [11], is only 0.005 parsecs. In order for protons from a cosmic-ray accelerator to reach us, the intervening magnetic field must connect us to the source and be coherent out to 100 parsecs since there are likely no sources within this distance. However, the direction of both regions is nearly perpendicular to the expected Galactic magnetic field direction [11] With nonstandard cosmic-ray diffusion, it is conceivable to account for these regions with a nearby cosmic-ray accelerator [12].

Another possibility is that one or both of the excesses could be caused by the heliosphere. This explanation is supported by the coincidence of region A with the direction of the heliotail ( RA74°, δ17° [13]), which is the direction opposite the motion of the solar system with respect to the local interstellar matter. The possibility that we are seeing neutron production in the gravitationally-focused tail of interstellar medium material has been considered and discarded in [12] because of insufficient target material.

In summary, Milagro has observed two unexplained regions of excess with high significance. Potential systematic causes have been examined and excluded. Both excesses are inconsistent with pure gamma rays with high confidence, and their energy spectra are moderately to strongly inconsistent with the spectrum of the isotropic cosmic-ray flux. In particular, the excess in region A can be modeled as hard spectrum protons with a cutoff.


  1. K. Nagashima, K. Fujimoto, and R. M. Jacklyn, J. Geophys. Res. 103, 17 429 (1998).
  2. D. L. Hall et al., J. Geophys. Res. 104, 6737 (1999).
  3. M. Amenomori et al., Science 314, 439 (2006).
  4. M. Ambrosio et al., Phys. Rev. D 67, 042002 (2003).
  5. K. Munakata et al., Phys. Rev. D 56, 23 (1997).
  6. B. E. Kolterman et al., 30th ICRC, Merida, Mexico (2007).
  7. R. Atkins et al. (Milagro Collaboration), Astrophys. J. 595, 803 (2003).
  8. T. P. Li and Y. Q. Ma, Astrophys. J. 272, 317 (1983).
  9. A. A. Abdo et al., Astrophys. J. 658, L33 (2007).
  10. M. Amenomori et al., in Turbulence and Nonlinear Processes in Astrophysical Plasmas, AIP Conf. Proc. No. 932 (AIP, New York, 2007), p. 283.
  11. J. L. Han, R. N. Manchester, A. G. Lyne, G. J. Qiao, and W. van Straten, Astrophys. J. 642, 868 (2006).
  12. L. Drury and F. Aharonian, Astropart. Phys. 29, 420 (2008).
  13. M. Witte, M. Banaszkiewicz, H. Rosenbauer, and D. McMullin, Adv. Space Res. 34, 61 (2004).

About the Authors

Image of A. A. Abdo
Image of B. Allen
Image of T. Aune
Image of D. Berley
Image of E. Blaufuss
Image of S. Casanova
Image of C. Chen
Image of B. L. Dingus
Image of R. W. Ellsworth
Image of L. Fleysher
Image of R. Fleysher
Image of M. M. Gonzalez
Image of J. A. Goodman
Image of C. M. Hoffman
Image of P. H. Hüntemeyer
Image of B. E. Kolterman
Image of C. P. Lansdell
Image of J. T. Linnemann
Image of J. E. McEnery
Image of A. I. Mincer
Image of P. Nemethy
Image of D. Noyes
Image of J. Pretz
Image of J. M. Ryan
Image of P. M. Saz Parkinson
Image of A. Shoup
Image of G. Sinnis
Image of A. J. Smith
Image of G. W. Sullivan
Image of V. Vasileiou
Image of G. P. Walker
Image of D. A. Williams
Image of G. B. Yodh

Related Articles

Synopsis: Helicons in a Lab Plasma
Plasma Physics

Synopsis: Helicons in a Lab Plasma

The plasma waves known as helicons can be created and measured in the laboratory even without confining walls. Read More »

Synopsis: Spatial Tests of Dark Matter

Synopsis: Spatial Tests of Dark Matter

Maps of merging galaxy clusters could help find signatures of dark matter based on its decay into photons. Read More »

Synopsis: Water Flow Helps Cells Move
Biological Physics

Synopsis: Water Flow Helps Cells Move

Water flowing through a cell’s membrane is essential to the process of changing cellular shape. Read More »

More Articles