Observation of Anisotropy of TeV Cosmic Rays with Two Years of HAWC

The study of the anisotropy in the arrival direction of cosmic rays has entered an era of precision measurement. Combined with improved modeling of the local interstellar medium, these measurements are maturing into an important way to simultaneously understand cosmic-ray acceleration and propagation. For a recent review, see Ahlers & Mertsch (2017).

Anisotropy is a well-studied consequence of standard propagation models where cosmic rays diffuse due to random magnetic fields and their sources are distributed inhomogeneously (Erlykin & Wolfendale 2006; Blasi & Amato 2012; Pohl & Eichler 2013; Sveshnikova et al. 2013). Anisotropy can also arise from motion relative to the rest frame of the cosmic rays (Compton & Getting 1935). Both scenarios result in a dominantly dipolar anisotropy, yet the predicted dipole amplitude is at least an order of magnitude larger than the observed value (Kumar & Eichler 2014; Mertsch & Funk 2015), and the measured dipole orientation cannot be explained by these simple models. Recent studies have included the effects of regular magnetic fields to probe the origins of the dipole direction (Ahlers et al. 2016) in the hopes of identifying the locations of accelerators contributing most to the locally observed cosmic-ray flux.

While the observed TeV cosmic-ray anisotropy is primarily dipolar with amplitude ∼10⁻³, it also contains smaller scale structure O (≲45°) with strength ∼10⁻⁴. It is likely that an initial dipolar distribution is distorted as it passes through the interstellar medium. For example, isotropic magnetic turbulence can create a feed down of angular power to higher multipole moments (Giacinti & Sigl 2012; Ahlers 2014; Ahlers & Mertsch 2015; Giacinti & Kirk 2017), and pitch-angle scattering (Giacinti & Kirk 2017) alters the shape of the large-scale multipoles. It is possible that heliospheric effects perturb the anisotropy, manifesting in the observed small-scale structure (Desiati & Lazarian 2013; Schwadron et al. 2014). Anisotropy can also result from nonstandard diffusion such as strong regularities in local magnetic field lines (Drury 2013; Harding et al. 2016). Thus, divergence from a pure dipole anisotropy provides a probe into the bulk properties of the interstellar medium.

Measuring anisotropy signals of O(10⁻³–10⁻⁴) at significant levels requires several key detector attributes. A large instantaneous sky coverage and long, uninterrupted observation periods are needed to achieve statistical uncertainties below the signal strength and to resolve features with large angular extent (180°). Only earthbound air-shower detectors fit these requirements, combining large fields of view, effective areas of ∼10⁴ m², and high duty cycles with long-term stability.

Cosmic-ray anisotropy has been observed in the energy range ∼500 GeV–100 PeV by air shower arrays such as Tibet-ASγ (Amenomori et al. 2005b, 2006, 2017), Milagro (Abdo et al. 2008, 2009), EAS-TOP (Aglietta et al. 2009), IceCube/IceTop (Abbasi et al. 2010, 2011, 2012; Aartsen et al. 2013, 2016), ARGO-YBJ (Bartoli et al. 2013, 2015; Di Sciascio 2013), and HAWC (Abeysekara et al. 2014). Below these energies, cosmic rays begin to follow geomagnetic field lines and no longer probe interstellar scales. Only the Pierre Auger Observatory has a significant measurement (Aab et al. 2017a, 2017b) above EeV energies. The most recent experimental overviews are given in Di Sciascio & Iuppa (2014) and in Ahlers & Mertsch (2017).

Air shower arrays must use the observations themselves to determine intrinsic detector acceptances, which limits sensitivity to the anisotropy component along the direction of the Earth's rotation, i.e., along the R.A. α in equatorial coordinates. These detectors are thus unable to recover a dipolar signal aligned with the equatorial poles, bounding the maximally recoverable dipole strength according to its orientation in decl. Direct modeling of the detector acceptance, as done in Aab et al. (2017a), can eliminate this bias but has not been demonstrated for arrays operating at TeV energies, as it requires an agreement between simulation and data over the full zenith angle range at the level of 10⁻⁵ or better, which is currently not achieved.

In this paper, we describe the results of an analysis of the cosmic-ray anisotropy on all angular scales and as a function of energy using the first two years of data recorded by the HAWC experiment. In a previous paper (Abeysekara et al. 2014), we used 113 days of data to study the small-scale anisotropy of cosmic rays above 1.7 TeV. Here, we extend on this analysis by also studying the large-scale anisotropy and using an improved energy estimator (Alfaro et al. 2017) to study the energy-dependence of both the small- and large-scale structures.

In addition, we apply a new iterative method (Ahlers et al. 2016) to reconstruct the maximally recoverable strength of the anisotropy. This method compensates for the reduction in the measured dipole strength caused by the fact that midlatitude detectors only see a fraction of the cosmic-ray dipole at any given time.

With the first two years of data taking, the HAWC array can currently study the cosmic-ray anisotropy up to energies of about 70 TeV. In future studies, we will extend the energy range to higher energies, using the same methods described in this paper.

This paper is organized as follows: we first describe the HAWC detector in Section 2, then the event selection and data set used for the measurement in Section 3. An explanation of the analysis methods is provided in Section 4. The results of the observed anisotropy are presented and discussed in Section 5. Section 6 summarizes the main conclusions of this work.

The High-Altitude Water Cherenkov (HAWC) Gamma-ray Observatory is an extensive air-shower array located at 4100 m a.s.l. on the slopes of Volcan Sierra Negra at 19°N in the state of Puebla, Mexico. While HAWC is designed to study the sky in gamma-rays between 500 GeV and 100 TeV, it is also sensitive to showers from primary cosmic rays up to multi-PeV energies.

The detector consists of a 22,000 m² array of 300 close-packed water Cherenkov detectors (WCDs), each containing 200 kiloliters of purified water and four upward-facing photomultiplier tubes (PMTs). As secondary air shower particles pass through the WCDs, the Cherenkov light produced is collected by the PMTs, permitting the reconstruction of primary particle properties including the local arrival direction, core location, and the energy. Further details on the HAWC detector can be found in Abeysekara et al. (2017).

The light-tight nature of the WCDs allows the detector to operate at nearly 100% up-time efficiency, with the data acquisition system recording air showers at a rate of ∼25 kHz. With a resulting daily sky coverage of 8.4 sr, HAWC is an instrument well-suited for measuring the cosmic-ray arrival direction distribution.

The data set for this study consists of 508 uninterrupted sidereal days between 2015 and 2017 May 1. The detector operated with 294 WCDs, recording about 1.2 × 10¹² air shower triggers. To determine the energy of the primary air shower particle, we apply a maximum-likelihood-based estimator that uses the lateral distribution of measured PMT signals as a function of simulated primary proton energy (Alfaro et al. 2017). To improve the estimated energy resolution, poorly reconstructed showers are removed from the data set by application of moderate event selection. The selection criteria are:

• 1.  
  Air shower events must pass a minimum multiplicity threshold of ≥75 PMTs. This improves angle and energy reconstruction accuracy.
• 2.  
  At least 1 PMT within 40 m of the core position (N_r40 ≥ 1) must record a signal. This criterion selects air showers landing on or near the array, resulting in a core resolution of better than 15 m above 10 TeV.
• 3.  
  The zenith angle acceptance range is 0°–60°.
• 4.  
  The data set is composed of periods covering complete sidereal days, hence events from incomplete days are not included. This removes nonuniformities in sky exposure along R.A., reducing systematics in the estimation of the reference map (see Section 4.1).

Requiring a multiplicity of ≥75 PMTs reduces the trigger rate to 23%. The remaining selection criteria further reduce the number of events by 52%, leaving a total of 123 billion air shower events.

Using the selection criteria and the likelihood energy estimator, we achieve ∼30% improvement in energy resolution compared to the multiplicity energy-proxy method from previous HAWC results (Abeysekara et al. 2014), and a ∼50% improvement over ARGO-YBJ (Di Sciascio 2013) below 10 TeV. This also permits more energy bins, as well as an increase in the median energy of the highest-energy bin. The median energy and 68% containment for the eight analysis bins are listed in Table 1 and shown in Figure 1 along with the energy bins from comparable experimental results.

Zoom In Zoom Out Reset image size

Table 1.  Reported Median Energy (With 68% Central Containment Region) and Fit of Two-dimensional Dipole Anisotropy (Amplitude and Phase) for Each Independent Energy Bin

Energy	Events	Amplitude	Phase	a1,1	a1,−1	${\chi }^{2}/{N}_{\mathrm{dof}}$
(TeV)		[×10−4]		[×10−4]	[×10−4]
$2.0({}_{+5.2}^{-1.4})$	4.0 × 1010	8.1 ± 0.4	429 ± 25	−17.1 ± 1.0	−15.9 ± 1.0	13.34
$3.0({}_{+7.1}^{-2.1})$	2.9 × 1010	8.9 ± 0.3	522 ± 20	−15.9 ± 0.9	−20.5 ± 0.9	1.24
$4.4({}_{+10.6}^{-3.2})$	2.4 × 1010	8.3 ± 0.3	456 ± 18	−16.7 ± 0.7	−17.1 ± 0.7	0.79
$6.8({}_{+14.0}^{-5.0})$	1.6 × 1010	10.1 ± 0.3	395 ± 16	−22.7 ± 0.8	−18.7 ± 0.8	0.85
$11.2({}_{+18.8}^{-7.9})$	7.9 × 109	11.9 ± 0.4	413 ± 19	−25.9 ± 1.1	−22.7 ± 1.1	0.81
$18.6({}_{+25.6}^{-12.7})$	3.8 × 109	13.8 ± 0.6	445 ± 24	−28.4 ± 1.6	−27.9 ± 1.6	1.07
$30.3({}_{+34.8}^{-19.3})$	1.8 × 109	14.4 ± 0.8	360 ± 32	−33.7 ± 2.3	−24.5 ± 2.3	1.25
$72.8({}_{+106.7}^{-44.9})$	1.6 × 109	6.7 ± 0.9	319 ± 73	−16.4 ± 2.5	−10.2 ± 2.5	1.03

Download table as:  ASCIITypeset image

The estimated energy exhibits a slight dependence on decl. as determined by simulation, shown in Figure 2. This is attributed to an increasing energy threshold with increasing zenith angle, as air showers must traverse more atmospheric overburden. The values presented in Figure 1 and Table 1 were calculated for the overhead sky (δ = 19°), thus this decl.dependent energy shift must be considered when viewing the resulting sky maps.

Zoom In Zoom Out Reset image size

Furthermore, for all maps in the lowest energy bin (2.0 TeV), there is a decreased range in decl. used for analysis, corresponding to the first two zenith bins (θ < 352) of the energy estimation described in Alfaro et al. (2017). This is due to the limited number of selected events having both large zenith angles and low reconstructed energies available to measure the anisotropy via the methods described in Section 4.

The estimated angular resolution given the selection criteria improves from 08 to 05 between 1 and 10 TeV, and plateaus at 05 above 10 TeV. A full description of the in situ angular resolution and energy-scale verifications with the cosmic-ray Moon shadow is presented in Alfaro et al. (2017).

The measured anisotropy requires comparison of the observed data event distribution ${ \mathcal D }$ with a background distribution or "reference map" ${ \mathcal B }$, which represents the detector response to an isotropic flux of cosmic rays. This reference map is not in itself isotropic because of effects of the detector exposure and geometry. In principle, the reference map can be obtained from a complete simulation of the detector response to an isotropic flux of cosmic rays, but as previously mentioned, measuring anisotropy at the 10⁻⁴ level requires an accuracy of the detector simulation that can currently not be achieved for detectors like HAWC. The reference map is therefore estimated from the data themselves.

We report the anisotropy distribution using two-dimensional sky maps, represented via the equal-area pixellation scheme provided by the HEALPix (Gorski et al. 2005) package. Maps are tessellated with 12 base pixels, and each pixel is further partitioned into N_side subdivisions. We chose a fine pixellation of N_side = 256, corresponding to a pixel width of 023 and solid angle of 1.6 × 10⁻⁵ sr.

4.1. Reference Map (Isotropic Expectation)

The reference map is determined using the method described in Atkins et al. (2003), due to its minimal assumptions on the data compared to similar background estimation methods (Amenomori et al. 2005a; Aglietta et al. 2009). This technique uses the number of events in local angular coordinates (the local detector acceptance) and in sidereal time (all-sky rate) to construct the expected counts map from an isotropic flux expectation.

For the calculation of the relative intensity of the cosmic-ray anisotropy, this paper uses for the first time a new analysis technique developed by Ahlers et al. (2016) that mitigates a common artifact of previous methods. For midlatitude detectors like HAWC, methods like that of Atkins et al. (2003) severely underestimate the relative intensity of any large-scale structure, in particular, the strength of the dipole component. The reason for this underestimation is the fact that these detectors have an instantaneous field of view that is much smaller than the size of the large-scale anisotropy structure. As a consequence, these large-scale structures are attenuated.

Over the course of a sidereal day, as the Earth rotates, the detector eventually accumulates an event distribution that shows the entire large-scale structure. However, a consequence of observing different parts of the anisotropy at different times is that the observed event distribution is a function of both the instantaneous detector exposure as well as the cosmic-ray anisotropy itself.

The new method overcomes the effect of the limited instantaneous exposure by simultaneously fitting for the anisotropy and the detector exposure, using a maximum-likelihood technique. The resulting equations cannot be solved explicitly, but best-fit solutions can be attained by iteratively adjusting the local acceptance and all-sky rate. The convergence criterion depends on the likelihood value of the calculated ${ \mathcal B }$ provided ${ \mathcal D }$ compared to the previous iteration. Typically, the calculation requires less than 20 iterations.

We demonstrate the stable convergence of the process with simulated dipoles of various orientations on the sky using a set of 10¹⁰ simulated events drawn from a HAWC-like sky exposure. The rate as a function of sidereal time was varied by a simple sinusoid of amplitude 5%. On top of the simulated events, dipoles of strength 10⁻³ with 10 different orientations in decl. (δ₀, from 0° to 90°) were added, providing 10 fake data sets. For each simulated data set, the differential relative intensity map was created using the iterative method, and the maximally recoverable dipole amplitude was then obtained. Figure 3 shows the stable convergence of the fit results after 1, 5, 10, and 20 iterations. We also simulated various sky coverages confirming that the fit obtains the maximally recoverable dipole amplitude.

Zoom In Zoom Out Reset image size

The all-sky rate varies by ∼5% due to diurnal pressure cycles in the upper atmosphere, and the local detector acceptance is verified to be stable for each sidereal day by evaluating the χ²-difference of local angular distributions in 2 minute intervals compared to the mean calculated over the entire sidereal day. With such minimal variations over the data set, we chose to sum the local detector acceptance and all-sky rate for all 508 sidereal days prior to the background calculation. Since the all-sky rate is binned in sidereal time bins of 1°, much smaller than the large- and small-scale features (>10°), variations within a single sidereal time bin have negligible impact on the observed features.

4.2. Relative Intensity and Significance Map

The amplitude of the measured anisotropy and its statistical strength are given by the sky maps in differential relative intensity δI and significance S, respectively. For a given pixel, δI is the fractional difference between the observed counts in that pixel, ${ \mathcal D }$, and the expected counts from the reference map, ${ \mathcal B }$:

$$\begin{eqnarray}&&\delta I=\displaystyle \frac{{ \mathcal D }}{{ \mathcal B }}-1.\end{eqnarray} \tag{ 1 }$$

The significance of δI is conservatively estimated via

$$\begin{eqnarray}&&S\simeq \sqrt{\displaystyle \frac{{ \mathcal D }}{1+{\alpha }_{\exp }}}\delta I,\end{eqnarray} \tag{ 2 }$$

based on Li & Ma (1983), where α_exp is the relative exposure of the data map compared to the reference map. The reference map is overexposed compared to the data because it uses information from all local pixels to calculate its values. The value of α_exp is found analytically via the method in Atkins et al. (2003), but a direct calculation was not determined for the iterative method. A conservative value of α_exp = 1 underestimates the statistical significance by at most 70%.

The results of the HAWC analysis are shown in Figure 4, which depicts the relative intensity of the arrival direction distribution of cosmic rays in equatorial coordinates for eight independent energy bins ranging from 2.0 to 72.8 TeV. Figure 5 shows the significance of the deviation from isotropy for the same energy bins, where negative values of significance correspond to pixels with δI < 1. The maps show significant deviation of the cosmic-ray flux from isotropy, dominated by a dipolar feature, which increases in strength up to 30.3 TeV while maintaining a nearly constant phase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Multipole Fitting

To better quantify the observed large-scale features, the relative intensity map δI is fit to the following truncated series of spherical harmonics:

$$\begin{eqnarray}&&\delta I=\displaystyle \sum _{{\ell }=1}^{{{\ell }}_{\max }}\displaystyle \sum _{m=-{\ell }}^{{\ell }}{a}_{{\ell }m}{Y}_{{\ell }m}(\theta ,\phi )\,\,(m\ne 0),\end{eqnarray} \tag{ 3 }$$

where ℓ_max is chosen to distinguish between large-scale and small-scale features. The choice of ℓ_max will be discussed in Section 5.

Since the local detector acceptance is estimated from the data, features in the anisotropy that only depend on decl. cannot be disentangled from decl.-dependent asymmetries in the detector acceptance. This results from the fact that the HAWC detector only samples the sky in the Earth's rotational direction. Thus, we set a_ℓ,0 = 0 for all ℓ, i.e., the m = 0 terms, to reflect the knowledge that the analysis method is insensitive to anisotropy oriented solely along the decl. direction.

In the fit, we chose to use the real-valued (tesseral) spherical harmonics,

$$\begin{eqnarray}{Y}_{{\ell }m}(\theta ,\phi )=\left\{\begin{array}{ll}\sqrt{\displaystyle \frac{2(2{\ell }+1)}{4\pi }\displaystyle \frac{({\ell }-| m| )!}{({\ell }+| m| )!}}{P}_{{\ell }}^{| m| }(x)\sin | m| \varphi & \,{\rm{if}}\ m\lt 0\\ \sqrt{\displaystyle \frac{2(2{\ell }+1)}{4\pi }\displaystyle \frac{({\ell }-m)!}{({\ell }+m)!}}{P}_{{\ell }}^{m}(x)\cos m\varphi & {\rm{if}}\ m\gt 0\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where $x=\cos \theta$, again discarding the m = 0 terms that are symmetric in R.A. The a_ℓm are then determined by a χ²-minimization fit of δI to Equation (3), and the variance for each pixel is calculated via propagation of uncertainties of the quantities ${ \mathcal D }$ and ${ \mathcal B }$ which comprise δI.

The strongest feature of the measured anisotropy is the dipole (l = 1), which can be more conveniently expressed as an amplitude and phase by projection onto R.A. The amplitude ${\widetilde{A}}_{1}$ can be expressed as the sum of the a_1,m terms added in quadrature:

$$\begin{eqnarray}&&{\widetilde{A}}_{1}=\sqrt{\displaystyle \frac{3}{8\pi }\displaystyle \sum _{m=1,-1}| {a}_{1m}{| }^{2}}\end{eqnarray} \tag{ 5 }$$

with variance

$$\begin{eqnarray}&&{\sigma }_{({\widetilde{A}}_{1})}^{2}=\displaystyle \frac{3}{8\pi {\sum }_{m=1,-1}| {a}_{1m}{| }^{2}}\displaystyle \sum _{m=1,-1}(| {a}_{1m}{| }^{2}{\sigma }_{({a}_{1m})}^{2}).\end{eqnarray} \tag{ 6 }$$

The maximally recoverable dipole amplitude ${\widetilde{A}}_{1}$ obtained via the iterative method is related to the true amplitude A₁ through the original decl. position of the maximum δ₀,

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \frac{{\widetilde{A}}_{1}}{\cos {\delta }_{0}}.\end{eqnarray} \tag{ 7 }$$

The term a₁₀ = 0 leaves one term that scales with the cosine (a_1,1) and one that scales with the sine (a_1,−1) of the dipole phase $\widetilde{{\phi }_{1}}$. This constrains the maximum amplitude to a decl. of 0°, simplifying the measured phase to

$$\begin{eqnarray}&&\widetilde{{\phi }_{1}}={\tan }^{-1}\left(\displaystyle \frac{{a}_{1,-1}}{{a}_{\mathrm{1,1}}}\right),\end{eqnarray} \tag{ 8 }$$

with variance

$$\begin{eqnarray}&&{\sigma }_{(\widetilde{{\phi }_{1}})}^{2}=\displaystyle \frac{1}{{\left({\sum }_{m=1,-1}| {a}_{1m}{| }^{2}\right)}^{2}}\displaystyle \sum _{m=1,-1}(| {a}_{1m}{| }^{2}{\sigma }_{({a}_{1m})}^{2}).\end{eqnarray} \tag{ 9 }$$

4.4. Angular Power Spectrum

The cosmic-ray anisotropy is not a pure dipole, and its full angular power spectrum reveals the strength of correlations at various angular scales. Figure 6 shows for each energy bin the pseudo-angular power spectra as derived from the anafast routine in HEALPix, where ${\tilde{C}}_{{\ell }}=\tfrac{1}{2{\ell }+1}{\sum }_{m=-{\ell }}^{{\ell }}| {a}_{{\ell }m}{| }^{2}$. The infinite series of multipoles was truncated at ℓ_max = 40 for quicker computation. Truncation at the maximum multipole that can be calculated for a HEALPix grid (ℓ_max = 3 N_side − 1) was done for comparison, with no noticeable effect on the reported power spectra.

Zoom In Zoom Out Reset image size

The uncertainties shown in Figure 6 are systematic, representing the 68% central containment region of the measured power spectrum for maps generated to have the same power spectrum as the data. The statistical uncertainties are smaller than the data points. These were determined by the 68% containment of angular power spectra derived from random sky maps drawn from a Poisson distribution using the true map as the mean.

The gray band shows the 90% confidence interval for the expected angular power spectrum for an isotropic map containing the same event statistics. The isotropic power band is flat across the multipole moments, and its magnitude level is determined by the number of events in the map. This band was derived in the same way as the statistical uncertainties, except the data fluctuated by a Poisson distribution was compared to the true data instead of the reference map. The deviation of the spectral points from this gray band represents the measured signal strength compared to isotropy.

4.5. Systematics

The resulting relative intensity maps, significance maps, and power spectra for each of the eight energy bins from 2.0 to 72.8 TeV are shown in Figures 4–6, respectively. The fit dipole amplitudes and phases obtained from these large-scale maps are shown in Figure 7. In order to enhance regional correlations, the sky maps have been smoothed by a circular top hat function of radius 10°, in accordance with other studies (Abdo et al. 2008; Abeysekara et al. 2014).

Zoom In Zoom Out Reset image size

Each map in Figure 4 shares the common significant features of having a broad region of deficit around α = 150° to α = 240° and a broad but more sharply peaked excess around α = 30° to α = 90°. The deficit grows in intensity with energy until the final bin at 72.8 TeV, where the feature diminishes. The center of the excess starts low in the HAWC field of view near δ = −15° and rises to about δ = −5° by 4.4 TeV, where it remains for the remaining energy bins. Its strength increases until 11.2 TeV, slowly diminishing until nearly disappearing in the 72.8 TeV bin. Starting at 6.8 TeV the excess develops an extension higher in decl. and slightly higher in R.A. This extension is strongest in the 30.3 TeV bin before also diminishing by 72.8 TeV. As shown in Figure 6, the dipole moments possess the most angular power for each bin, also reflected in the evolution of the broad deficit. The second strongest moment is the quadropole, whose power remains fairly constant at around 1.5 × 10⁻⁷ for all energies save the highest-energy bin. The octupole moment is typically half of the quadrupole moment. In most bins there is a rapid decrease in power from these first three moments to the sextupole (ℓ = 4), at which point we differentiate between the large and small angular scales. Significant anisotropy is seen up to ℓ = 10 (characteristic angular scale of 18°) until statistics dip below 5 billion events above 18.6 TeV. The decrease in ${\tilde{C}}_{{\ell }}$ as a function of ℓ, especially for ℓ > 3 becomes more rapid with increasing energy, and the anisotropy becomes less significant due to the rising noise floor of maps with fewer data. This noise level is represented as the expected power spectrum of an isotropic cosmic-ray distribution, shown by the gray bands in Figure 6. The large uncertainty for the ℓ = 1 term in the lowest energy bin results from the reduced decl. range available for the multipole fit. As the integrated sky coverage decreases, angular power from lower multipoles becomes increasingly degenerate with power from higher multipoles.

A summary of the dipole fit parameters obtained per the methods of Section 4.3, and the median cosmic-ray energies and numbers of events for each bin are given in Table 1. The dipole component is detected at a significant level in all bins, and as shown in Figure 7, its amplitude steadily increases with energy from 8.1 × 10⁻⁴ to 14.4 × 10⁻⁴ until the final bin at 72.8 TeV where its value drops to 6.7 × 10⁻⁴. As shown in Figure 5, this bin also has the most significant excess near α = 300°, the decl. of the Cygnus region, which has more than one extended TeV gamma-ray emission feature. Determining the precise contribution of gamma-ray contamination will be considered in future studies.

The resulting phases and amplitudes from the dipole fits are compared with measurements from other experiments in Figure 8, with the HAWC measurements (green squares) being in fair agreement with the observed trends. Though not all phases are consistent within statistical errors, systematics that are not accounted for in the estimation of the reference maps may contribute to the 5°−10° differences between HAWC and ARGO-YBJ (black diamonds) at energies below 10 TeV. The phases measured by HAWC and IceCube are consistent within the overlapping energy range, noting that the reference map methods used here and in the IceCube study are nearly identical.

Zoom In Zoom Out Reset image size

The evolution with energy of the fit amplitudes matches well with previous results, showing a steady rise until a sudden decrease between 50 and 100 TeV. For the HAWC measurement, the highest-energy bin at 72.8 TeV has nearly the same number of events as at 30.3 TeV, yet its dipole amplitude is reduced by more than half. The energy scale of the amplitude behavior is in slight disagreement with several other experiments. This tension could be resolved by shifting along the abscissa, as the experiments's energy scales may be offset relative to one another. For example, the energy scale reported for HAWC in this work may vary by 5% (Alfaro et al. 2017), while the ARGO-YBJ proton energy is reported to within 13% (Bartoli et al. 2015).

IceCube also shows a slight discrepancy in the amplitude scale of the anisotropy with other measurements, including HAWC. It is possible that differences in chemical compositions at the detector level are the cause, as according to simulations, IceCube measures a higher-rigidity composition than IceTop (Aartsen et al. 2016) and potentially other ground-based air shower detectors. While higher-rigidity particles follow field lines more closely, an anisotropy signal from many parsecs away might be distorted and diminished by the nearby magnetic fields of the Earth and the Sun. For reference, using the simulated composition described in Alfaro et al. (2017), we find that the fraction of events passing the selection criteria from proton and helium primaries decreases with energy from 91% at 3 TeV to 78% at 10 TeV, and to 69% at 100 TeV.

In Figure 8, only the HAWC measurement shown by the green squares uses a two-dimensional fit as described in Section 4.3. All previous measurements represent fits of the dipole component to a Fourier series after projection of the relative intensity sky map onto a single decl. band. To match the method presented by the other experimental results, we also include one-dimensional fit dipole parameters (purple squares), being between 75% and 85% of the two-dimensional fit values. This suggests that previous dipole amplitudes also are underestimated, primarily affecting results from detectors with larger integrated fields of view (e.g., HAWC, ARGO-YBJ, and Tibet) as compared to those with more limited fields of view such as IceCube.

In addition to the large-scale structure, the sky maps have significant angular power at small angular scales (ℓ ≥ 4) as shown in the angular power spectra, and in the relative intensity and significance maps after subtraction of the ℓ ≤ 3 fit multipoles, shown in Figures 9 and 10. The three most significant regions of excess previously observed with Milagro (Abdo et al. 2008) and HAWC (Abeysekara et al. 2014) are Regions A, B, and C (referred to as Regions 1, 2, and 4 by ARGO-YBJ Bartoli et al. 2013). These are more apparent in the relative intensity and significance maps of Figure 11, where all events have been combined into a single map having median energy of $2.6({}_{+9.9}^{-2.0})$ TeV to study each region's morphology. The excess defined as Region 3 by ARGO-YBJ (Bartoli et al. 2013) has a maximum significance of 4.6σ at α = 2512, δ = 195 in the combined HAWC map, with a mean relative intensity for the entire region of (1.4 ± 0.3) × 10⁻⁴, approximately 14% greater than that measured by ARGO-YBJ. Figure 12 shows zoomed-in views of each regional excess observed by HAWC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Regions A and C are characterized as relatively symmetric excesses with an extent of ≈30°, as compared to Region B, which has a similar width in R.A. but is elongated by nearly a factor of two in decl. Similar morphology was observed for all three regions in the previous small-scale study by HAWC (Abeysekara et al. 2014). A significant feature previously not observed by HAWC nor ARGO-YBJ is labeled Region D in Figure 11. This new excess, also shown in Figure 13, occupies about 25° in decl. while also being elongated by about a factor of two in R.A., though is considerably weaker than the other regions, having a maximum relative intensity of (1.7 ± 0.2) × 10⁻⁴.

Region A is the most prominent feature in δI as seen in Figure 9, being present at all energies, while the shape of Region B is evident up to 18.6 TeV and that of Regions C and D are difficult to discern given the color scale. The presence of these features in significance as a function of energy is shown in Figure 10, with Regions A and B being significant up to 18.6 TeV, Region D peaking at 4.4 TeV and Region C being strongest below 4.4 TeV.

Table 2 presents for the four regions the locations in (α, δ) of the peak significances (σ_max) and the corresponding δI values found in the all-bin combined maps. The relative intensity values are consistent with those found for Regions A, B, and C at σ_max in Abeysekara et al. (2014), while the locations of σ_max are shifted by ∼5° from their previously identified values. The relative intensity spectra as a function of energy extracted from the maps in Figure 9 at these locations are shown in Figure 14, including that of Region D. Regions B and C are characterized by relatively flat spectra across the energy bins as compared to that of Region D, which increases in δI up to the 4.4 TeV bin before decreasing, and Region A, which has a positive slope up to 6.8 TeV. This is consistent with the previous spectral measurement of Region A by HAWC (Abeysekara et al. 2014) as shown in the right panel of Figure 14. Due to improved uncertainties from a data set an order of magnitude greater, the shape of the spectrum from this study is much more constrained, especially for E ≥ 10 TeV. As was done in Abeysekara et al. (2014), we estimate the statistical significance of the spectral slope of Region A by comparing a linear fit of δI(log E) to similar fits performed across the field of view. The distribution of slopes across the sky (excluding points within 20° of Regions A, B, C, and D) follows a Gaussian distribution with a mean of zero and width 0.69 × 10⁻⁴. The slope fit at the location of Region A is (1.66 ± 0.46) × 10⁻⁴, falling 2.4σ from the all-sky mean.

Zoom In Zoom Out Reset image size

From Figures 4 and 9 it appears that the large-scale structure as well as Region A are approximately oriented along the local interstellar magnetic field ${{\boldsymbol{B}}}_{\mathrm{LIMF}}$ inferred from Interstellar Boundary Explorer (IBEX) measurements. In equatorial coordinates (α, δ), this direction is (23443 ± 069, 163 ± 045), and (485 ± 069, −212 ± 045) for $-{{\boldsymbol{B}}}_{\mathrm{LIMF}}$ (Zirnstein et al. 2016). This alignment is consistent with local conditions playing a role in shaping the observed cosmic-ray anisotropy, providing insight into the structure of the local interstellar medium and the heliospheric environment (Desiati & Lazarian 2013; Schwadron et al. 2014).

The HAWC Observatory has observed significant cosmic-ray anisotropy on both large and small scales using 123 × 10⁹ events comprising one of the largest TeV anisotropy data sets to date. Implementing an energy estimation technique that has been verified using the cosmic-ray Moon shadow (Alfaro et al. 2017), we have achieved an unprecendented energy resolution for measuring the energy-dependence of the anisotropy over eight analysis bins. Additionally, a new maximum-likelihood method was used to recover a minimally biased estimate of the expected intensity of an isotropic signal.

The energy dependence of the large-scale phase and amplitude is found to be consistent with observations made by other detectors in the northern hemisphere. Similarly, the morphology and relative intensity spectra of the three most significant small-scale regions of excess are consistent with previous measurements made by HAWC (Abeysekara et al. 2014) and ARGO-YBJ (Bartoli et al. 2013, 2015; Di Sciascio 2013). Finally, due to the increased statistics of the current data set, it is possible to further constrain the spectrum of Region A from previous Milagro and HAWC measurements.

Along with continued optimization of the event selection, the ever-growing HAWC data set will facilitate increasingly accurate descriptions of the anisotropy as a function of energy, providing additional insights into the nature of local accelerators and the interstellar environment. Furthermore, the novel techniques used in this study allow for collaboration with other observatories using data sets consisting of targeted cosmic-ray energy bands. Combination of HAWC data with the IceCube cosmic-ray data set is an ongoing effort that will form a nearly complete map of the cosmic-ray sky at TeV energies.

Dedicated to the memory of our honorable colleague, beloved mentor, and dear friend Stefan Westerhoff. We acknowledge support from the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México (grants 271051, 232656, 260378, 179588, 239762, 254964, 271737, 258865, 243290, 132197, 281653), Laboratorio Nacional HAWC de rayos gamma; L'OREAL Fellowship for Women in Science 2014; Belgian American Educational Foundational Fellowship; Wallonie-Bruxelles International; Red HAWC, México; DGAPA-UNAM (grants IG100317, IN111315, IN111716-3, IA102715, 109916, IA102917); VIEP- BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant DEC-2014/13/B/ST9/945; Coordinación de la Investigacíon Científica de la Universidad Michoacana. We thank Markus Ahlers for help with the implementation of the maximum-likelihood method. We gratefully acknowledge Scott DeLay for his dedicated efforts in the construction and maintenance of the HAWC experiment. Additional thanks to Luciano Díaz and Eduardo Murrieta for technical support.

Software: HEALPix (Gorski et al. 2005).