# Discovery of Localized Regions of Excess 10-TeV Cosmic Rays

Phys. Rev. Lett. 101, 221101

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 $60 m×80 m$ pond surrounded by a sparse $200 m×200 m$ 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 $∼2 sr$ field of view, operates with a $>90%$ duty cycle, and has a trigger rate from cosmic rays of $∼1700 Hz$, 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].

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 $∼1 TeV$ 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].

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$.

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° and $10°<δ<20°$. Region $B$ is defined as the union of two boxes: $117° and $15°<δ<40°$, and $131° 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 $log⁡e(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 $log⁡10(Ec/GeV)=4.0-0.5(stat)+0.4$ for region $A$ and $log⁡10(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 $log⁡e(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 ( $RA≈74°$, $δ≈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.

## References

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).

## Related Articles

Cosmology

### Viewpoint: Cosmic Clues from Mini Clumps of Dark Matter

Searches for ultracompact clumps of cold dark matter have come up empty, but these nondetections place new limits on the early expansion history of the Universe. Read More »

Superfluidity

### Synopsis: Doubling Up with Superfluids

Researchers mixed two superfluids of different atoms together and observed that vortices in one affected those in the other—evidence of mutual interaction between the two species. Read More »

Materials Science

### Viewpoint: Improving Electronic Structure Calculations

A new approach to calculating the properties of molecules and solids may offer higher accuracy at reasonable computational cost, accelerating the discovery of useful materials. Read More »