All-Sky Measurement of the Anisotropy of Cosmic Rays at 10 TeV and Mapping of the Local Interstellar Magnetic Field

TL;DR

This study provides the first full-sky measurement of cosmic ray anisotropy at 10 TeV, revealing the local interstellar magnetic field's direction and properties through combined data from HAWC and IceCube.

Contribution

It introduces a comprehensive all-sky analysis of cosmic ray anisotropy at 10 TeV, combining data from two observatories to reduce bias and map the local interstellar magnetic field.

Findings

Determined the horizontal dipole components of anisotropy.

Mapped the interstellar magnetic field direction.

Estimated the vertical dipole component of anisotropy.

Abstract

We present the first full-sky analysis of the cosmic ray arrival direction distribution with data collected by the HAWC and IceCube observatories in the Northern and Southern hemispheres at the same median primary particle energy of 10 TeV. The combined sky map and angular power spectrum largely eliminate biases that result from partial sky coverage and holds a key to probe into the propagation properties of TeV cosmic rays through our local interstellar medium and the interaction between the interstellar and heliospheric magnetic fields. From the map we determine the horizontal dipole components of the anisotropy $\delta_{0h} = 9.16 \times 10^{-4}$ and $\delta_{6h} = 7.25 \times 10^{-4}~(\pm0.04 \times 10^{-4})$. In addition, we infer the direction ($229.2\pm 3.5^\circ$ RA , $11.4\pm 3.0^\circ$ Dec.) of the interstellar magnetic field from the boundary between large scale excess and deficit regions from which we estimate the missing corresponding vertical dipole component of the large scale anisotropy to be $\delta_N \sim -3.97 ^{+1.0}_{-2.0} \times 10^{-4}$.

Tables

Table 1. Table 1: Comparison of the IceCube and HAWC datasets. The median primary particle energy, angular resolution and number of remaining events is shown in the sub-columns after applying only quality cuts and after applying both energy and quality cuts. The angular resolution of IceCube corresponds to the DST dataset that relies on faster, less accurate reconstructions as well as less stringent quality cuts. The energy cuts applied are chosen to lower the median energy of IceCube data from 20 TeV down to 10 TeV. In the case of HAWC, the cuts are aimed to raise the median energy of HAWC data from 2 TeV up to 10 TeV.

	IceCube		HAWC
Latitude	90∘ S		19∘ N
Detection method	muons produced by CR		air showers produced by CR and $\gamma$
Field of view	-90∘/-16∘ ($\delta$), $\sim$4 sr (same sky over 24h)		-30∘/68∘ ($\delta$), $\sim$2 sr (8 sr observed/24 h)
Livetime	1742 days over a period of 1826 days		519 days over a period of 653 days
Detector trigger rate	2.5 kHz		25 kHz
	Quality cuts	Energy and quality cuts	Quality cuts	Energy and quality cuts
Median primary energy	$20$ TeV	$10$ TeV	$2$ TeV	$10$ TeV
Approx. angular resolution	$2^{\circ }-3^{\circ }$	$2^{\circ }-6^{\circ }$	${0.4}^{\circ }-0^{\circ }.8$	${0.4}^{\circ }-1^{\circ }.0$
Events	$2.8\times {10}^{11}$	$1.7\times {10}^{11}$	$7.1\times {10}^{10}$	$2.8\times {10}^{10}$

Table 2. Table 2: Relative mass composition for 10 TeV median energy cosmic-rays in the two samples as determined from CORSIKA Monte Carlo simulations  (Heck et al., 1998 ) weighted to a Polygonato spectrum  (Jörg R. Hörandel, 2003 ) . Errors reflect statistical uncertainties in the simulation datasets.

	IceCube (10 TeV)	HAWC (10 TeV)
Proton	0.756 $\pm$ 0.018	0.6160 $\pm$ 0.0054
He	0.195 $\pm$ 0.009	0.3110 $\pm$ 0.0014
CNO	0.028 $\pm$ 0.004	0.0467 $\pm$ 0.0004
NeMgSi	0.013 $\pm$ 0.002	0.0191 $\pm$ 0.0001
Fe	0.008 $\pm$ 0.002	0.0078 $\pm$ 0.0001

Table 3. Table 3: Spherical Harmonic Coefficients [ 10 − 4 superscript 10 4 10^{-4} ]

$\ell$	$m=$1	2	3	4	5	6	7
1	-13.26 +10.49$i$
2	-0.21 -3.6$i$	-7.20 +2.05$i$
3	1.75 -1.7$i$	-2.03 +0.13$i$	0.20 +0.17$i$
4	1.70 -0.52$i$	0.07 +1.69$i$	-0.86 -0.8$i$	-1.19 +0.04$i$
5	0.58 +0.27$i$	-0.07 -1.1$i$	-1.64 -0.051$i$	0.18 -0.15$i$	-0.11 -1.5$i$
6	0.80 -0.88$i$	-0.24 -0.38$i$	-0.10 +0.63$i$	0.13 -1.2$i$	0.27 +0.47$i$	1.65 -0.53$i$
7	0.44 -0.67$i$	0.37 +0.15$i$	-0.21 -0.14$i$	-0.70 +0.04$i$	0.84 -0.27$i$	0.13 -0.54$i$	0.07 +0.91$i$
8	0.26 +0.06$i$	0.14 -0.47$i$	-0.39 -0.22$i$	-0.42 +0.72$i$	-0.15 -0.15$i$	-0.72 -0.61$i$	0.42 +0.36$i$
9	0.11 -0.88$i$	-0.29 -1.3$i$	0.22 -0.17$i$	0.12 -0.56$i$	-0.01 -0.34$i$	0.60 +0.47$i$	-0.06 -0.48$i$
10	0.21 -0.97$i$	0.25 -0.5$i$	0.21 -0.65$i$	0.09 -0.088$i$	-0.10 +0.12$i$	0.11 -0.017$i$	0.02 +0.19$i$
11	0.56 -0.39$i$	0.06 -0.42$i$	-0.15 -0.68$i$	-0.04 +0.05$i$	-0.26 +0.04$i$	-0.07 -0.26$i$	-0.16 +0.25$i$
12	0.40 +0.07$i$	0.19 -0.56$i$	-0.27 -0.48$i$	-0.17 -0.1$i$	-0.13 -0.18$i$	-0.03 -0.23$i$	0.33 +0.13$i$
13	0.45 -0.33$i$	-0.04 -0.69$i$	0.17 -0.92$i$	-0.26 -0.6$i$	0.13 +0.24$i$	-0.08 +0.02$i$	0.04 +0.04$i$
14	0.57 -0.16$i$	0.13 -0.53$i$	0.17 -1.1$i$	-0.31 -0.089$i$	0.08 -0.09$i$	-0.25 -0.12$i$	-0.05 +0.22$i$
$\ell$	$m=$8	9	10	11	12	13	14
8	-0.54 +0.19$i$
9	0.15 +0.64$i$	-0.04 +0.45$i$
10	0.22 +0.12$i$	-0.66 -0.57$i$	-0.26 +0.38$i$
11	0.25 +0.02$i$	-0.21 -0.4$i$	0.15 -0.25$i$	-0.06 -0.18$i$
12	0.37 +0.09$i$	-0.46 +0.25$i$	-0.13 +0.20$i$	-0.08 +0.21$i$	0.04 -0.18$i$
13	0.11 +0.13$i$	0.13 -0.13$i$	-0.35 -0.098$i$	0.39 +0.45$i$	-0.01 -0.3$i$	0.41 -0.17$i$
14	-0.13 +0.34$i$	0.36 -0.11$i$	-0.04 -0.072$i$	-0.11 -0.17$i$	-0.19 +0.32$i$	0.13 +0.21$i$	0.18 +0.35$i$

Table 4. Table 4: Magnetic field alignment. The last two rows correspond to measurements of the large-scale anisotropy from this study. The R.A. measurement in the last row is obtained from the dipole vector and the decl. is obtained from the ℓ = 2 ℓ 2 \ell=2 quadrupole component. The second to last column corresponds to the angular distance Δ ​ ψ Δ 𝜓 \Delta\psi between the boundary fit and the various LIMF estimates. The last column gives the corresponding vertical dipole component under the assumption that the dipole is oriented towards the given decl. Error in parentheses for dipole and quadruple correspond to systematic uncertainties.

Source	R.A. [∘]	decl. [∘]	$\Delta \psi$ [∘]	${\delta }_N$ [${10}^{-4}$]
Funsten et al. (2013)	$229.7\pm 1.9$	$23.3\pm 2.5$	9.0	-5.03
Frisch et al. (2015)	$237.9\pm 16$	$22.2\pm 16$	12.2	-4.77
Zirnstein et al. (2016)	$234.4\pm 0.7$	$16.3\pm 0.6$	6.5	-3.42
Boundary Fit	$229.2\pm 3.5$	$11.4\pm 3.0$	–	-2.36
Dipole/Quadrupole	$218.4\pm 0.3$ ($\pm 2.6$)	$20.7\pm 0.3$ ($\pm 2.6$)	–	-4.41

Equations

$$\mathcal{P}=\frac{1}{N}\sum _{i=0}^N\frac{({\zeta }_i-⟨{\zeta }_i⟩{)}^2}{{{\sigma }_{{\zeta }_i}}^2}$$

$$\delta I_{\mathfrak{a}}=\frac{N_{\mathfrak{a}}-⟨N_{\mathfrak{a}}⟩}{⟨N_{\mathfrak{a}}⟩}.$$

$$S_i=\sqrt{2}{\left(-{\mu }_{i,\mathrm{on}}+{\mu }_{i,\mathrm{off}}+n_i\log \frac{{\mu }_{i,\mathrm{on}}}{{\mu }_{i,\mathrm{off}}}\right)}^{1/2}.$$

$$n_i$$

$${\mu }_{i,\mathrm{on}}$$

$${\mu }_{i,\mathrm{off}}$$

$$\delta I({\mathbf{u}}_i)=\sum _{\ell =1}^{\infty }\sum _{m=-\ell }^{\ell }{a}_{\ell m}{Y}_{\ell m}({\mathbf{u}}_i).$$

$$\mathbf{\delta }\equiv ({\delta }_{0h},{\delta }_{6h},{\delta }_N)=\sqrt{\frac{3}{2\pi }}(-\Re (a_{11}),\Im (a_{11}),a_{10}),$$

$$({\delta }_{0h},{\delta }_{6h})=({\tilde{A}}_1\cos {\alpha }_1,{\tilde{A}}_1\sin {\alpha }_1),$$

$${\mathcal{C}}_{\ell }=\frac{1}{2\ell +1}\sum _{m=-\ell }^{\ell }|a_{\ell m}{|}^2,$$

$${\tilde{\mathcal{C}}}_{\ell }=\frac{1}{2\ell }\sum _{m=-\ell ,m\ne 0}^{\ell }|a_{\ell m}{|}^2,$$

$$\frac{\Delta I}{I}=(\gamma +2)\frac{v}{c}\cos ({\theta }_v),$$

$$\Delta \{\ln R\}=\beta \cdot \Delta P$$

$$\mathcal{E}(t,\phi ,\theta )\simeq E(t)A(\phi ,\theta ),$$

$${\mu }_{\tau i}\simeq I_{\tau i}{\mathcal{N}}_{\tau }{\mathcal{A}}_i,$$

$$\mathcal{L}(n|I,\mathcal{N},\mathcal{A})=\prod _{\tau i}\frac{({\mu }_{\tau i}{)}^{n_{\tau i}}{e}^{-{\mu }_{\tau i}}}{n_{\tau i}!},$$

$$\lambda =\frac{\mathcal{L}(n|I,\mathcal{N},\mathcal{A})}{\mathcal{L}(n|I^{(0)},{\mathcal{N}}^{(0)},{\mathcal{A}}^{(0)})},$$

$${\mathcal{N}}_{\tau }^{(0)}$$

$${\mathcal{A}}_i^{(0)}$$

$$\sum _i{\mathcal{A}}_i=1.$$

$$\mathcal{E}(t,\phi ,\theta )\simeq \sum _{\mathrm{sector}s}{E}^s(t){\mathcal{A}}^s(\phi ,\theta ).$$

$$I_{\mathfrak{a}}^{\star }$$

$${\mathcal{N}}_{\tau }^{s\star }$$

$${\mathcal{A}}_i^{s\star }$$

Full text

astroparticle physics, cosmic rays, ISM: magnetic fields

††journal: ApJ††software: HEALPix/healpy (version 1.9.1, Górski et al. (2005)), CORSIKA (version 7.40, Heck et al. (1998)), ROOT (version 6.04/12, Brun & Rademakers (1996)), Matplotlib (version 1.5.0, Hunter (2007)), Astropy (version 1.1, Astropy Collaboration et al. (2013); The Astropy Collaboration et al. (2018)), SciPy (version 0.16.1, http://www.scipy.org/), NumPy (version 1.11.1,Oliphant (2015)), Python programming language (Python Software Foundation, https://www.python.org/), PolSpice (version 3.0.3, http://www2.iap.fr/users/hivon/software/PolSpice).

Email: juan.diazvelez@alumnos.udg.mx

Email: analysis@icecube.wisc.edu

††thanks: Earthquake Research Institute,

University of Tokyo, Bunkyo, Tokyo 113-0032, Japan

1 Introduction

A number of theoretical models predict an anisotropy in the distribution of arrival directions of cosmic rays that results from the distribution of sources in the Galaxy and diffusive propagation of these particles (Erlykin & Wolfendale, 2006; Blasi & Amato, 2012; Ptuskin, 2012; Pohl & Eichler, 2013a; Sveshnikova et al., 2013; Kumar & Eichler, 2014a; Mertsch & Funk, 2015). Although the observed distribution of cosmic rays is highly isotropic, several ground-based experiments located either in the northern or southern hemisphere have observed small but significant variations in the arrival direction distribution of TeV to PeV cosmic rays with high statistical accuracy, in both large and medium angular scales (Nagashima et al., 1998; Hall et al., 1999; Amenomori et al., 2005, 2006; Guillian et al., 2007; Abdo et al., 2008, 2009; Aglietta et al., 2009; Munakata et al., 2010; Abbasi et al., 2010, 2011; De Jong, 2011; Abbasi et al., 2012; Aartsen et al., 2013; Bartoli et al., 2013; Abeysekara et al., 2014; Bartoli et al., 2015; Aartsen et al., 2016; Amenomori et al., 2017; Bartoli et al., 2018; Abeysekara et al., 2018b). The observed large-scale anisotropy has an amplitude of about $10^{-3}$ and small-scale structures of amplitude of $10^{-4}$ with angular size of $10^{\circ}$ to $30^{\circ}$.

For previously reported measurements that rely on time-integrated methods (Alexandreas et al., 1993; Atkins et al., 2003), a difference between the instantaneous and integrated field of view of the experiments can lead to an attenuation of structures with angular size larger than the instantaneous field of view (Ahlers et al., 2016). For this analysis, we apply an optimal reconstruction method that can recover the amplitude of the projected large-scale anisotropy. The limited integrated field of view of the sky in all of these individual measurements also makes it difficult to correctly characterize such an anisotropy in terms of its spherical harmonic components and produce a quantitative measurement of the large scale characteristics, such as its dipole or quadrupole component, without a high degree of degeneracy (Sommers, 2001). The resulting correlations between the multipole spherical harmonic terms $a_{\ell m}$ bias the interpretation of the cosmic ray distributions in the context of particle diffusion in the local interstellar medium (LISM). In this joint analysis by the High-Altitude Water Cherenkov (HAWC) and IceCube collaborations we have combined data from both experiments at 10 TeV median primary particle energy to study the full-sky anisotropy. Important information can be obtained from the power spectrum of the spherical harmonic components at low $\ell$ (large scale), which is most affected by partial sky coverage. It should be noted that neither observatory is sensitive to variations across decl. bands since events recorded from a fixed direction in the local coordinate system can only probe the cosmic-ray flux in a fixed decl. band $\delta$. As a result, the dipole anisotropy can only be observed as a projection onto the equatorial plane. However, some information about the vertical component can be inferred from medium- and small-scale structures.

2 The HAWC Detector

The HAWC Gamma-Ray Observatory is an extensive air-shower detector array located at 4100 m a.s.l. on the slopes of Volcán Sierra Negra at $19^{\circ}$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 with an instantaneous field of view of about $2$ sr.

The detector consists of a 22,000 m2 array of 300 close-packed water Cherenkov detectors (WCDs), each containing 200 metric tons of purified water and four upward-facing photomultiplier tubes (PMTs). At the bottom of each WCD, three 8-inch Hamamatsu R5912 photomultiplier tubes (PMTs) are anchored in an equilateral triangle of side length 3.2 meters, with one 10-inch high-quantum efficiency Hamamatsu R7081 PMT anchored at the center.

As secondary air shower particles pass through the WCDs, the produced Cherenkov light 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, 2018a).

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 $\sim$ 25 kHz. With a resulting daily sky coverage of $8.4$ sr and an angular resolution of $0^{\circ}.4$ for energies above 10 TeV, HAWC is an ideal instrument for measuring the cosmic-ray arrival direction distribution with unprecedented precision.

3 The IceCube Detector

The IceCube Neutrino Observatory, located at the geographic South Pole, is composed of a neutrino detector in the deep ice and a surface air-shower array. The in-ice IceCube detector consists of 86 vertical strings containing a total of 5,160 optical sensors, called digital optical modules (DOMs), frozen in the ice at depths from 1,450 meters to 2,450 meters below the surface of the ice. A DOM consists of a pressure-protective glass sphere that houses a 10-inch Hamamatsu photomultiplier tube together with electronic boards used for detection, digitization, and readout. The strings are separated by an average distance of 125 m, each one hosting 60 DOMs equally spaced over the kilometer of instrumented length. The DOMs detect Cherenkov radiation produced by relativistic particles passing through the ice, including muons and muon bundles produced by cosmic-ray air showers in the atmosphere above IceCube. These atmospheric muons form a large background for neutrino analyses, but also provide an opportunity to use IceCube as a large cosmic-ray detector. Further details on the IceCube detector can be found in Aartsen et al. (2017).

All events that trigger IceCube are reconstructed using a likelihood-based method that accounts for light propagation in the ice (Ahrens et al., 2004). The fit provides a median angular resolution of $3^{\circ}$ according to simulation (Abbasi et al., 2011) but worsens past zenith angles of approximately $70^{\circ}$ Aartsen et al. (2017). This is not to be confused with the $\sim 0^{\circ}.6$ angular resolution of IceCube for neutrino-induced tracks of where more sophisticated reconstruction algorithms and more stringent quality cuts are applied. The energy threshold of cosmic-ray primaries producing atmospheric muons in IceCube is limited by the minimum muon energy required to penetrate the ice. As a result, the primary particle energy threshold increases with larger zenith angles as muons must travel increasingly longer distances through the ice. This is accounted for in the analysis as described in Sec. 4. Due to the limited data transfer rate available from the South Pole, cosmic-ray induced muon data are stored in a compact data storage and transfer (DST) format (Abbasi et al., 2011), containing the results of the angular reconstructions described as well as some limited information per event. However, detailed information such as PMT waveforms used for these reconstructions is not kept. The preliminary reconstructions encoded in the DST rely on faster, less accurate methods than those applied to the filtered dataset used in most neutrino analyses.

4 The Dataset

The dataset selected for this analysis is composed of 5 years of data collected by the IceCube Neutrino Observatory between May 2011 and May 2016, as well as 2 years of data from the HAWC Gamma-Ray Observatory collected between May 2015 and May 2017. In order to reduce bias from uneven exposure along R.A., only full sidereal days of continuous data-taking were chosen for this study. The residual contribution of the dipole anisotropy induced by the motion of the Earth around the Sun is estimated to be on the order of $10^{-5}$, which is smaller than the statistical error of this analysis (see section 7.2). Cuts are applied to each dataset in order to improve the angular resolution and energy resolution of reconstructed events. In the case of HAWC these include a cut on the number of active optical sensors in order to increase the information available for the reconstruction of the shower. A cut on the reconstructed zenith angle excludes events with $\theta>57^{\circ}$ where the quality of reconstructions decreases rapidly. A cut is also applied on the variable CxPE40 which corresponds to the effective charge measured in the PMT with the largest effective charge at a distance of more than 40 m from the shower core with CxPE40 $>40$. The effective charge $Q_{\mathrm{eff}}$ scales the charge of higher-efficiency central 10-inch PMTs by a factor of 0.46 relative to the 8-inch PMTs so that all optical sensors are treated equally. The value of CxPE40 is typically large for a hadronic events (Abeysekara et al., 2017). Finally, in order to identify and exclude gamma-ray candidates, a cut is applied on $\mathcal{P}$, that describes the “clumpiness” of the air shower (Abeysekara et al., 2017) with $\mathcal{P}>1.8$. $\mathcal{P}$ is defined using the lateral distribution function of the air shower. $\mathcal{P}$ is computed using the logarithm of the effective charge $\zeta_{i}={\rm{log}}_{10}(Q_{\rm{eff},i})$. For each PMT hit, $i$, an expectation is assigned $\langle\zeta_{i}\rangle$ by averaging the $\zeta_{i}$ in all PMTs contained in an annulus containing the hit, with a width of 5 meters, centered at the core of the air shower. $\mathcal{P}$ is then calculated using the $\chi^{2}$ formula:

[TABLE]

The errors $\sigma_{\zeta_{i}}$ are assigned from a study of a sample strong gamma-ray candidates in the vicinity of the Crab nebula. The $\mathcal{P}$ variable essentially requires axial smoothness.

In the case of IceCube we apply a cut on the reduced likelihood of the directional reconstruction ($\mathrm{RlogL}<15$), defined as the best-fit log-likelihood divided by the number of degrees of freedom in the fit (Ahrens et al., 2004) which gives an estimate of the goodness of fit for the angular reconstruction. There is also a cut on the number of direct photoelectrons and the corresponding length of the track $N_{\mathrm{dir}}>9\cos(\theta)$ and $l_{\mathrm{dir}}>200\cos(\theta)$ meters. This cut depends on the reconstructed zenith angle $\theta$ in order to preserve sufficient statistics near the horizon. Photons are considered direct when the time residual (i.e., the delay in their arrival time due to scattering in the ice) falls within a time window of -15 ns to +75 ns with respect to the geometrically expected arrival time from the reconstructed track (Ahrens et al., 2004).

Table 1 shows the characteristics of both experiments next to each other. The two detectors have different energy responses and this results in a difference in the median energy. In order to select data that are consistent between the two detectors, we have applied additional cuts on the reconstructed energy of events: in the case of HAWC we use an energy reconstruction based on the likelihood method described in Alfaro et al. (2017) to select events with reconstructed energies at or above 10 TeV. In the case of IceCube we apply a cut in the two dimensional plane of number of hit optical sensors (which act as a proxy for muon energy) and the cosine of the reconstructed zenith angle, as described in Abbasi et al. (2012). As a result of the overburden of ice described in Sec. 3, for a given number of hit optical sensors, events at larger zenith angles are produced by cosmic-ray particles with higher energy (Abbasi et al., 2012; Aartsen et al., 2016). The energy resolution is primarily limited by the relatively large fluctuations in the fraction of the total shower energy that is transferred to the muon bundle and is of the order of 0.5 in $\log_{10}(E/\mathrm{GeV})$ (Aartsen et al., 2016).

Figure 1 shows the distribution of data as a function of decl. The resulting energy distribution of the two datasets is shown in Figure 2. As a result of the applied energy cuts, both cosmic-ray data sets have a median primary particle energy of approximately 10 TeV with little dependence on zenith angle (Figure 3). The energy response of the observatories covers a 68% range of approximately 3 TeV - 40 TeV, in the case of IceCube, and 2.5 TeV - 30 TeV for HAWC around the median energy.

The two experiments have different response to the cosmic ray mass composition. This is largely due to the detection method. Particles entering Earth’s atmosphere (15 to 20 km above sea level) interact with nuclei in air and produce a cascade of secondary particles. This particle cascade continues to grow until ionization becomes the dominant energy loss mechanism. The depth Xmax at which this happens depends on both the energy of the primary particle, and its mass. Lighter nuclei penetrate deeper than heavier nuclei. As a result, the altitude of extended air shower arrays such as HAWC can affect the response of the detector to different nuclei since they are sensitive to the electromagnetic component of the particle shower. In contrast, the IceCube in-ice detector observes cosmic rays through the detection of deep penetrating muons produced from the decay of charged pions and kaons generated in the early interactions. As a a result, for the same composition, IceCube’s response to different cosmic-ray nuclei differs from that of HAWC.

If the first interaction occurs at a lower air density (and higher elevation), mesons are more likely to decay to muons (and neutrinos) instead of re-interacting and producing lower energy pions and other secondary particles. As a result, the two experiments react differently to changes in atmospheric temperature and pressure.

The data from both experiments are dominated by light nuclei (protons and alpha particles) as can be seen in Table 2. All of the cuts applied were chosen based on CORSIKA Monte Carlo simulations (Heck et al., 1998) weighted to a Polygonato spectrum (Jörg R. Hörandel, 2003) and detailed simulations of the detector response.

5 Analysis

We compute the relative intensity as a function of J2000 equatorial coordinates ($\alpha$, $\delta$) by binning the sky into an equal-area grid with a bin size of $0^{\circ}.9$ using the HEALPix library (Górski et al., 2005). The angular distribution can be expressed as $\phi(\alpha,\delta)=\phi^{\mathrm{iso}}I(\alpha,\delta)$, where $\phi^{\mathrm{iso}}$ corresponds to the isotropic flux (i.e., the flux averaged over the full celestial sphere), and $I(\alpha,\delta)$ is the relative intensity of the flux as a function of R.A. $\alpha$ and decl. $\delta$ in celestial coordinates. Given that cosmic rays have been observed to be mainly isotropic, the flux is dominated by the isotropic term and therefore the anisotropy $\delta I=I-1$ is small.

The relative intensity gives the amplitude of deviations in the number of counts $N_{\mathfrak{a}}$ from the isotropic expectation $\langle N_{\mathfrak{a}}\rangle$ in each angular bin $\mathfrak{a}$. The residual anisotropy $\delta I$ of the distribution of arrival directions of the cosmic rays is calculated by subtracting a reference map that describes the detector response to an isotropic flux

[TABLE]

In order to produce this reference map, we must have a description of the arrival direction distribution if the cosmic rays arrived isotropically at Earth.

Ground-based experiments observe cosmic rays indirectly by detecting the secondary air shower particles produced by collisions of the cosmic-ray primary in the atmosphere. The observed large-scale anisotropy has an amplitude of about $10^{-3}$ but our simulations are not sufficiently accurate to describe the detector response at this level. We therefore calculate this expected flux from the data themselves in order to account for detector dependent rate variations in both time and viewing angle. For Earth-based observatories, such a method requires averaging along each decl. band, thus washing out the vertical dependency (i.e. as a function of decl. $\delta$) in the relative intensity map $\delta I_{\mathfrak{a}}$.

A common approach is to estimate the relative intensity and detector exposure simultaneously using time-integration methods (Alexandreas et al., 1993; Atkins et al., 2003). However, these methods can lead to an under- or overestimation of the isotropic reference level for detectors located at mid latitudes, since a fixed position on the celestial sphere is only observable over a relatively short period every day. As a result, the total number of cosmic ray events from this fixed position can only be compared against reference data observed during the same period. Therefore, time-integration methods can strongly attenuate large-scale structures exceeding the size of the instantaneous field of view (Ahlers et al., 2016).

5.1 Maximum Likelihood Method

For this analysis, we have relied on the likelihood-based reconstruction described in Ahlers et al. (2016) and recently applied in the study of the large-scale cosmic-ray anisotropy by HAWC (Abeysekara et al., 2018b). The method does not rely on detector simulations and provides an optimal anisotropy reconstruction and the recovery of the large-scale anisotropy projected on to the equatorial plane for ground-based cosmic ray observatories located in the middle latitudes as HAWC. The generalization of the maximum likelihood method for combined data sets from multiple observatories that have exposure to overlapping regions of the sky is described in Appendix A.

5.2 Statistical Significance

In order to calculate the statistical significance of anisotropy features in the final reconstructed map, Ahlers et al. (2016) generalizes the method in Li & Ma (1983) to account for the optimization process of the time-dependent exposure. The significance map (in units of Gaussian $\sigma$) is then calculated as

[TABLE]

For each pixel $i$ in the celestial sky, we define expected on-source and off-source event counts from neighbor pixels in a disc of radius $r$ centered on that pixel. For this analysis we have chosen a radius of 5∘. Given the set of pixels $\mathcal{D}_{i}$, the observed and expected counts are

[TABLE]

where $\mathcal{A}_{\tau j}$ is the relative acceptance of the detector in pixel $j$ and sidereal time bin $\tau$, $\mathcal{N}_{\tau}$ gives the expected number of isotropic events in sidereal time bin $\tau$, $I_{j}$ is the relative intensity, and where $I=I^{\rm reference}+I^{\rm residual}$ is divided into a contribution from the reference map and the residual relative intensity. For small-scale features, $I^{\rm reference}$ corresponds to the first 3 spherical harmonic components ($\ell\leq 3$) of the relative intensity. In order to distinguish excess and deficit, we multiply Eq. 3 by the sign of each smoothed pixel $\delta I_{i}$ in the anisotropy map.

5.3 Harmonic Analysis and Dipole Fit

The relative intensity can be decomposed as a sum over spherical harmonics $Y^{\ell m}$,

[TABLE]

The vector components of the dipole in terms of the spherical harmonic expansion $Y_{\ell m}$ in equatorial coordinates are related to the $a_{\ell m}$ coefficients with

[TABLE]

where $\Re(a_{11})$ and $\Im(a_{11})$ are respectively, the real and imaginary components of $a_{11}$, and taking into account that $a_{1-1}=-a^{*}_{11}$ and $a_{10}=a^{*}_{10}$ (see Ahlers & Mertsch (2017)).

From equation 8 and the $a_{\ell m}$ coefficients, one can obtain the horizontal components of the dipole $\delta_{0h}$ and $\delta_{6h}$ with respect to the [math]h and $6$h R.A. axes. The phase and amplitude of the projected dipole on the equatorial plane are given by

[TABLE]

where $\alpha_{1}$ is the phase and $\tilde{A}_{1}$ is the amplitude of the projected dipole on the equatorial plane and it is related to the true amplitude $A_{1}$ through the dipole inclination $\delta_{0}$ with $\tilde{A}_{1}=A_{1}\cos\delta_{0}$.

5.4 Angular Power Spectrum

The angular power spectrum for the relative intensity field is defined as:

[TABLE]

for each value of $\ell$. Since this analysis is not sensitive to the vertical component of the anisotropy, the largest recoverable dipole amplitude $\tilde{A}$ has the terms $m=0$ missing and we can only measure a pseudo power spectrum $\tilde{\mathcal{C}}_{\ell}$:

[TABLE]

The angular power spectrum provides an estimate of the significance of structures at different angular scales of $\sim$ 180${}^{\circ}/\ell$. In the ideal case of a $4\pi$ sky coverage, the multipole moments ${a}_{\ell m}$ of the reconstructed anisotropy would carry all the information of the anisotropy (except for the $m=0$ vertical component terms). However, as will be discussed in Section 7.2, partial sky coverage of individual experiments further limits the amount of information that can be obtained from the reconstructed pseudo multipole moment spectrum.

6 Results

The measured relative intensity map is shown in Figure 4. A smoothing procedure was applied to all maps using a top-hat function in which a single pixel’s value is the average of all pixels within a 5∘ radius. The map shows an anisotropy in the distribution of arrival directions of cosmic rays with 10 TeV median primary particle energy that extends across both hemispheres. The significance (smoothed by summing over pixels) of the IceCube region reflects the much larger statistics in the IceCube dataset compared to that from HAWC at energies of $\sim$10 TeV.

Figure 5 is the residual small-scale anisotropy after subtracting the fitted multipole from the spherical harmonic expansion with $\ell\leq 3$ from the large-scale map in Figure 4 in order to reveal structures smaller than $60^{\circ}$. The large-scale structure and significant small-scale structures in Figures 4 and 5 are largely consistent with previous individual measurements, as shown in Figure 6. Observed features extend across the horizon of both datasets. The one referred to as “region A” by the Milagro Collaboration (Abdo et al., 2008) roughly extends from $(54^{\circ},-16^{\circ})$, to $(78^{\circ},18^{\circ})$, in equatorial coordinates ($\delta$, $\alpha$). The so called “region B” (Abdo et al., 2008) corresponds to the boundary between the excess and deficit regions (see Figure 4) in the northern sky that appears as a small scale feature (see Figure 5) for short integration times.

We obtain the $a_{\ell m}$ through a transformation of spherical harmonics using the HEALPix function map2alm. The results are presented in Table 3. The horizontal components of the dipole obtained from equation (8) using the $a_{\ell m}$ values in Table 3 are $\delta_{0h}=9.16\times 10^{-4}$ and $\delta_{6h}=7.25\times 10^{-4}~{}(\pm 0.04\times 10^{-4})$, respectively, with respect to the [math]h and $6$h R.A. axes. The dipole amplitude and phase $\tilde{A}_{1}=(1.17\pm.01)\times 10^{-3}$, $\alpha_{1}=38.4\pm 0^{\circ}.3$ measured in this combined study are shown in Figure 6 along with previously published results from other experiments in the TeV-PeV primary particle energy range. The combined systematic uncertainty in the amplitude and phase of the dipole are expected to be $\delta\tilde{A}_{1}\sim 0.06\times 10^{-3}$, and $\delta\alpha_{1}\sim 2^{\circ}.6$ respectively (see section 7).

The angular power spectrum for the combined dataset in Figure 7 provides an estimate of the significance of structures at different angular scales of $\sim$ 180${}^{\circ}/\ell$. Biases are substantially reduced with the likelihood method and by eliminating degeneracy between multipole moments with a nearly full sky coverage. The angular power spectrum can therefore be considered to be the physics fingerprint of the observed 10 TeV anisotropy, providing information about the propagation of cosmic rays and the turbulent nature of the Local Interstellar Magnetic Field (LIMF) (Giacinti & Sigl, 2012; Ahlers & Mertsch, 2017). The large discrepancy between the combined and individual datasets is the result of the limited sky coverage by each experiment. This systematic effect will be discussed in Section 7.2. A residual limitation in this analysis is the fact that ground-based experiments are generally not sensitive to the vertical component of the anisotropy as discussed by Abeysekara et al. (2018b) and Ahlers et al. (2016), as mentioned earlier.

The measured quadrupole component has an amplitude of $6.8\times 10^{-4}$ and is inclined at $20.7~{}\pm 0^{\circ}.3$ above (and below) equatorial plane. As with the dipole, the fitted quadrupole component from the spherical harmonic expansion is also missing the $m=0$ terms. However, the combination of $a_{21}$ and $a_{22}$ non-vertical quadrupole components can still provide valuable information. The experimental determination of the vertical components of the anisotropy would require accuracies better than the amplitude of the anisotropy ($\sim 10^{-3}$). This becomes easier at ultra-high energies where a dipole of much larger amplitude has been observed (Aab et al., 2017). The full-sky coverage also provides better constraints for fitting the $\ell=2$ and $\ell=3$ multipole components and reduces correlations between spherical harmonic expansion coefficients $a_{\ell m}$.

7 Systematics Studies

7.1 Overlapping Region

We have studied two adjacent $\delta$ bands at -20∘ for HAWC and IceCube data near the horizon of each detector (see Figure 8). The HAWC band extends from -21∘ to -19∘ while the IceCube band extends from -22∘ to -20∘. The large structure between the two datasets is consistent though small structures differ. It is worth noting that the overlap region is where we expect to find the largest difference in median energy between the two datasets (see Figure 3). The angular resolution of both detectors also decreases toward the horizon. While HAWC data has a smaller point spread function at this decl. and is sensitive to structures on smaller scales, IceCube has better statistics so the structures are more significant. One particular feature that stands out is the excess in HAWC around $\alpha=50^{\circ}$ that coincides with the so called “region A”. There appears to be a corresponding small excess in the IceCube data. It is also worth noting that statistics in this region are quickly decreasing with increasing zenith angle as is the quality of angular reconstructions. As a result, $\delta$ bins closer to the horizon contain a high level of contamination from bins in higher zenith angles.

7.2 Partial Sky Coverage

Incomplete coverage of the sky leads to an underestimation of the angular power of the dipole perpendicular to the axis of rotation of the Earth. The pseudo-moments of the projected dipole, $a_{11}$ and $a_{1-1}$, are corrected by a geometric factor introduced by Ahlers et al. (2016) in order to estimate the true moments $\hat{a}_{11}$ and $\hat{a}_{1-1}$. Furthermore, there is a degeneracy between different $\ell$ pseudo-modes under partial sky coverage that primarily affects the multipolar components $\ell=2$, $\ell=3$, and to a lesser degree, $\ell=4$ as has been previously studied by Sommers (2001). This effect is evident in Figure 9 which corresponds to a dipole injected horizontally in the direction $\delta_{\mathrm{6h}}$. The partial coverage of the sky produces an artificial quadrupole, octupole and hexadecapole that, in the case of a horizontal dipole, decrease in power with greater celestial coverage. The horizontal axis indicates the maximum observable decl. $\delta_{\mathrm{max}}$, keeping $\delta_{\mathrm{min}}=-90^{\circ}$.

From Figure 9 it is possible to see that the spurious quadrupole and octupole components (which are significant for partial integrated sky coverage) are reduced to an amplitude to order $10^{-5}$ in this analysis. Figure 10 shows the correlation matrix (Efstathiou, 2004) of the different $\ell$-modes up to $\ell=30$ calculated using the PolSpice111 PolSpice website: http://www2.iap.fr/users/hivon/software/PolSpice/. software package. The correlation between $\ell$-modes due to partial sky coverage is appreciable for larger $\ell$, though to a lesser degree.

7.3 Seasonal Variations and Local Variations in Solar Time

The relative motion of the Earth around the Sun can introduce a systematic solar dipole, a dipole anisotropy analogous to the Compton-Getting effect (Compton & Getting, 1935) produced by the motion of Earth around the Sun, that points in the direction of Earth’s orbital velocity vector. The influence of diurnal variations (such as the solar dipole) on the sidereal anisotropy can be estimated from the influence it has on the anti-sidereal distribution in a frame with 364.24 cycles per year (see, e.g. Guillian et al. (2007)). Any significant variations in this frame result from a modulation of the solar frame and represents a systematic effect of the solar frame on the sidereal anisotropy (Aartsen et al., 2016). The anti-sidereal distribution of the HAWC dataset has a maximum amplitude of $5\times 10^{-5}$. Both contamination from the solar dipole and atmospheric pressure variations are included in this systematic. For IceCube, the same systematic uncertainty is at the level of $\sim 3\times 10^{-5}$. The worst-case uncertainty on the reconstructed phase of the dipole is $\delta\alpha=2^{\circ}.6$ and a combined systematic uncertainty of $\delta\tilde{A}=6\times 10^{-5}$ for the dipole amplitude.

The solar dipole anisotropy produced by the motion of Earth around the Sun is given by the equation

[TABLE]

where $I$ is the cosmic-ray intensity, $\gamma$ is the index of the differential energy spectrum of cosmic rays, $v$ is the velocity of the Earth, $c$ is the speed of light and $\theta_{v}$ is the angle between the direction of the reconstructed cosmic rays and the direction of the velocity vector (Compton & Getting, 1935). This vector rotates by $360^{\circ}$ such that, after one year, the effect is ideally completely cancelled for 100% duty cycle of observation. However, a residual dipole can be introduced if the data does not cover an integer number of years with uniform coverage. In other words, any gaps in data taking can result in a slight bias to the measured dipole. A solar dipole anisotropy at the level of $10^{-4}$ has been previously observed at several TeVs (Amenomori et al., 2004, 2006; Abdo et al., 2009; Abbasi et al., 2011, 2012; Bartoli et al., 2015). Based on Monte Carlo studies, the residual contribution solar dipole that results from gaps in data taking is estimated to be of order $\sim 10^{-5}$ for the HAWC dataset, which is smaller than the statistical error of this analysis. In the case of IceCube, the detector has an uptime of 99% (see  Aartsen et al. (2017)) reduced to an uptime of 95.4% after selecting full sidereal days. As a result, the systematic effect of data gaps is smaller (Abbasi et al., 2012).

In addition to variations caused by the anisotropy and the solar dipole, there may also be local variations in the detection of cosmic rays caused by changes in atmospheric conditions, such as pressure and temperature, and also by changes in the detector. For 10 TeV energies, HAWC is located below the shower maximum Xmax for all primary masses. As a result, an increase in pressure leads to an increase of the atmospheric overburden which results in an attenuation of shower sizes. Atmospheric overburden is related to ground pressure $p$ as $X_{0}=p/g$, where $g=$ 9.87 m s*-2* is the local gravitational acceleration (Abbasi et al., 2013). In first order approximation, the simple correlation between the change in the logarithm of the rate $\Delta\{\ln R\}$ and the surface pressure change $\Delta P$ is

[TABLE]

where $\beta$ is the barometric coefficient (Tilav et al., 2010). The variations in atmospheric pressure at the HAWC site are primarily due to atmospheric tides driven by temperature and a small contribution from gravitational tides (Zhang et al., 2010). We have studied the effect of atmospheric pressure variations by applying a correction to the data rate to account for measured changes in pressure at the HAWC site. The procedure involves determining the correlation coefficient between the surface pressure data and the detector rate from Eq. 13 in order to weight individual events. This yields a barometric coefficient of $\beta=-0.0071$ hPa*-1* The residual contamination from atmospheric variations is estimated to be on the order of $\sim 10^{-6}$. Temperature variations in the stratosphere can introduce a similar effect with a 24h cycle and a 365 day cycle. However, this effect is small for latitudes near the equator and in the case of the daily variations, it is a smaller effect than that of pressure variations.

In contrast with HAWC, where the event rate is anti-correlated with atmospheric pressure and with the effective temperature of the stratosphere, the muon rate in IceCube is directly correlated with the effective temperature (Tilav et al., 2010). Event rate variations in IceCube have an annual period since one day at the South Pole lasts 365 days instead of 24 hours. In the case of IceCube there are also faster atmospheric variations of lower amplitude but these approximately affect the event rate globally in all azimuth directions (with a maximum Kolmogorov-Smirnov distance below $9.6\times 10^{-4}$ for daily variations at a 90% confidence level). Due to the geometry of the detector and its location at the South Pole, this also means that such variations equally affect every angle in R.A.

Seasonal variations in the effective temperature can introduce modulations in the intensity of the Solar dipole. As a result, the Solar dipole would not average to zero over a full year and thus would produce a residual bias. However, the amplitude of the anti-sidereal distribution indicates that this is not a significant effect.

8 Discussion

The combined sky map of arrival direction distribution of the 10 TeV cosmic rays collected by HAWC and IceCube and the corresponding power spectrum of its spherical harmonics components, may provide important hints on the origin of the observation. In particular, the angular power spectrum can reveal information about how cosmic rays propagate through the interstellar medium while the large-scale arrival direction distribution provides hints about the structure of the nearby LIMF and the heliosphere.

8.1 Cosmic Ray Propagation in the Interstellar Medium

The angular power spectrum in Figure 7 shows two different regimes: a steeply falling slope at large scales $\ell=1,2,3$ and a softer slope at small scales $\ell>3$. This suggests that two different mechanisms are responsible for the structures observed in the sky map. The steep portion of the angular power spectrum may be associated with large scale diffusive processes (over many mean free paths) across the interstellar medium, as suggested by Erlykin & Wolfendale (2006); Blasi & Amato (2012); Ptuskin (2012); Pohl & Eichler (2013b); Sveshnikova et al. (2013); Savchenko et al. (2015); Ahlers (2016); Giacinti & Kirk (2017). On the other hand, the softer slope portion appears to be consistent with non-diffusive pitch angle scattering effects on magnetic turbulence within the mean free path (Giacinti & Sigl, 2012) and with that obtained from numerical calculations of sub-PeV protons propagating through incompressible magnetohydrodynamic turbulence (López-Barquero et al., 2016). In Ahlers (2014) it is shown that under certain conditions, those small-scale structures arise as natural consequence of hierarchical evolution of angular scales under Liouville’s theorem.

The dipole component of the anisotropy may provide a hint into the direction of the large scale cosmic ray density gradient on the equatorial plane, thus linking the observed anisotropy with possible contributions of the closest sources, such as the Vela supernova remnant at a distance of about 0.3 kpc and with an age of about 11 kyr (Ahlers & Mertsch, 2017). The fact that Vela is located within the large-scale excess region of the sky is consistent with it being a potential source contributing to the large-scale anisotropy. However, predictions of the anisotropy amplitude depend on many unknown factors such as the actual contributing source (or sources), the diffusion coefficient, and the unknown component of the anisotropy perpendicular to the equatorial plane that complicate such calculations.

The measured amplitude and phase in this study is consistent with observations from multiple experiments that show a turning point in the energy dependency of the dipole component amplitude at an energy scale of 10 TeV (see Figure 6). After initially increasing with energy, the dipole amplitude begins to decrease above 10 TeV, while the phase has an abrupt change at the 100 TeV energy scale where the amplitude begins to increase again. Cosmic rays with rigidity of 10 TV have a gyro-radius of about 700 AU in a 3 $\mu$G magnetic field, which is comparable to the transversal size of the heliosphere (i. e. perpendicular to the long axis) (Pogorelov et al., 2009, 2013; Pogorelov, 2016). It is reasonable to assume that at lower energies the heliospheric influence is important, while above 10 TV the interstellar influence is progressively more important (Desiati & Lazarian, 2013). An understanding of how interstellar propagation of 10 TV-scale cosmic rays influences the arrival direction distribution must, therefore, also take into account heliospheric effects (Schwadron et al., 2014; López-Barquero et al., 2017; Zhang & Pogorelov, 2016). An alternative approach is to study cosmic ray anisotropy above 100 TV rigidity (Aartsen et al., 2013, 2016), where the heliospheric influence is expected to be negligible. In this case, the arrival direction distribution can be used to probe the global properties of interstellar turbulence by fitting theoretical models to observations (Giacinti & Kirk, 2017). However, at high energies a full-sky study is currently not possible with the dataset used in this analysis due to limited statistics.

8.2 Large-scale Anisotropy and the Local Interstellar Magnetic Field

Figure 11 shows the direction of the LIMF from Zirnstein et al. (2016) and the corresponding equator (the continuous black line), the so-called $B-V$ plane defined by the LIMF and the direction of the Sun’s velocity through the interstellar medium $v_{ISM}$, as well as the direction of the velocity relative to the local standard of rest $v_{LSR}$. The figure also shows the location of the Geminga and Vela supernova remnants as possible contributing sources, and those of the Cygnus X-1 X-ray binary and Galactic center GC for reference. The location of the Galactic plane is shown as a red line. A fit to the plane defined by the small-scale feature that marks the boundary between the excess and deficit regions ($\sim 115^{\circ}$ R.A.) is shown in Figure 12. The fit yields a vector pointed towards $(\alpha_{\mathrm{fit}},\delta_{\mathrm{fit}})=(229.2\pm 3^{\circ}.5,11.4\pm 3^{\circ}.0)$ in J2000 equatorial coordinates, as shown in Figure 11, along with the corresponding equator (the crossed black curve). The direction is located $9^{\circ}$ from the LIMF inferred by the Interstellar Boundary Explorer (IBEX) from the emission of energetic neutral atoms (ENA) originating from the outer heliosphere (see  Funsten et al. (2013). This point is also located $6^{\circ}.5$ from the LIMF direction reported by Zirnstein et al. (2016) and consistent with the average LIMF direction obtained from the polarization of stars within 40 pc by Frisch et al. (2015). This is shown in Figure 13 and summarized in Table 4 along with the value of $\alpha$ obtained from the dipole fit and the value of $\delta$ obtained from the quadrupole fit. The errors on the fit are derived from the $\chi^{2}$ distribution shown in Figure 12 and don’t include possible systematics uncertainties from the missing $m=0$ dipole component.

The fact that the dipole component of the full-sky cosmic ray anisotropy map is approximately aligned with the direction of the LIMF (or at least its projection on the equatorial plane) is probably not a coincidence, since we expect diffusion to be anisotropic with the fastest propagation along the magnetic field lines (Effenberger, F. et al., 2012; Kumar & Eichler, 2014b; Schwadron et al., 2014; Mertsch & Funk, 2015). Assuming that the observed dipole points in this direction, it is possible to estimate the amplitude of the vertical component. The measured amplitude of the horizontal component of the dipole $\tilde{A}_{1}$ is related to the true amplitude $A_{1}$ through the dipole inclination $\delta_{0}$ with $\tilde{A}_{1}=A_{1}\cos\delta_{0}$, from which we obtain a value for the vertical dipole vector component of $\delta_{N}=\tilde{A}_{1}\tan{\delta_{0}}\sim-3.97^{+1.0}_{-2.0}\times 10^{-4}$ for the various magnetic field assumptions (see Table 4).

If we assume that the dipole component must be aligned with the LIMF, the observed deviation could be explained as due to the relative motion of the observer with respect to a frame in which the cosmic ray distribution is isotropic, called the Compton-Getting Effect (Compton & Getting, 1935; Gleeson & Axford, 1968). The heliosphere could also have a significant warping effect on 10 TeV cosmic ray arrival direction distribution, mostly due to the LIMF draping curvature around the heliosphere (Pogorelov et al., 2009). As a result, the dipole component of the cosmic ray anisotropy could be out of alignment from the LIMF. Future studies, with full-sky maps at different particle rigidities, could provide a more powerful tool to probe the properties of the interstellar and heliospheric magnetic fields.

9 Conclusions

We have used experimental data collected by the HAWC Gamma Ray Observatory and the IceCube Neutrino Observatory to compile, for the first time, a nearly full-sky map of the arrival direction distribution of cosmic rays with median energy of 10 TeV. The combined analysis accounts for the difference in instantaneous and time-integrated field of view of the HAWC observatory and provides an integrated field of view that extends from $-90^{\circ}$ to $+76^{\circ}$ in decl. The almost full-sky observation eliminates the degeneracy between the spherical harmonic components and provides a tool to probe the properties of particle diffusion in the interstellar medium and of interstellar magnetic turbulence. The corresponding angular power spectrum suggests that two different mechanisms are responsible for the observed angular scale features. The ordering of cosmic ray anisotropy along the LIMF is supported by fitting the boundary between deficit and excess, which points to the direction $(\alpha_{\mathrm{fit}},\delta_{\mathrm{fit}})=(229.2\pm 3^{\circ}.5,11.4\pm 3^{\circ}.0)$ that is consistent with various observations. We obtained the phase and amplitude of the dipole component projected onto the equatorial plane to be $\tilde{A}_{1}=(1.17\pm.01)\times 10^{-3}$, $\alpha_{1}=38.4\pm 0^{\circ}.3$. Based on the assumption that the true dipole is aligned along the LIMF, we estimated the missing vertical component to be $\delta_{N}\sim-3.97^{+1.0}_{-2.0}\times 10^{-4}$.

Acknowledgements

Dedicated to the memory of our honorable colleague, and dear friend Stefan Westerhoff. The IceCube collaboration acknowledges the significant contributions to this manuscript from M. Ahlers, P. Desiati., and J. C. Díaz-Vélez. The HAWC Collaboration acknowledges additional contributions from D. W. Fiorino.

The HAWC Collaboration acknowledges the 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)(Cátedras 873, 1563, 341), Laboratorio Nacional HAWC de rayos gamma; L’OREAL Fellowship for Women in Science 2014; Red HAWC, México; DGAPA-UNAM (grants IG100317, IN111315, IN111716-3, IA102715, IN109916, IA102917, IN112218); 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, DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society - Newton Advanced Fellowship 180385. Thanks to Scott Delay, Luciano Díaz and Eduardo Murrieta for technical support.

The IceCube Collaboration acknowledges the support from: USA – U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison, Open Science Grid (OSG), Extreme Science and Engineering Discovery Environment (XSEDE), U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle physics computational facility at Marquette University; Belgium – Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany – Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden – Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia – Australian Research Council; Canada – Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark – Villum Fonden, Danish National Research Foundation (DNRF), Carlsberg Foundation; New Zealand – Marsden Fund; Japan – Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea – National Research Foundation of Korea (NRF); Switzerland – Swiss National Science Foundation (SNSF).

Appendix A Generalized Maximum likelihood method for multiple observatories with overlapping fields of view

The method developed by Ahlers et al. (2016) assumes that the detector exposure $\mathcal{E}$ per solid angle and sidereal time $t$ accumulated over many sidereal days can be expressed as a product of its angular-integrated exposure $E$ per sidereal time and relative acceptance $A$ (normalized as $\int d\Omega A(\Omega)=1$):

[TABLE]

with the assumption that the relative acceptance of the detector does not strongly depend on sidereal time.

For each observatory, the number of cosmic rays expected from an angular element $\Delta\Omega_{i}$ of the local coordinate sphere corresponding to coordinates $(\theta_{i},\varphi_{i})$ in a sidereal time interval $\Delta t_{\tau}$ is

[TABLE]

where $\mathcal{N}_{\tau}\equiv\Delta t_{\tau}\phi^{\rm iso}{E}(t_{\tau})$ gives the expected number of isotropic events in sidereal time bin $\tau$ independent of pixel, $\mathcal{A}_{i}$ is the relative acceptance of the detector for pixel $i$, and $I_{\tau i}\equiv{I}({\bf R}(t_{\tau}){\bf n}^{\prime}(\Omega_{\mathfrak{i}}))$ is the relative intensity observed in local coordinates during time bin $\tau$. ${\bf R}(t){\bf n}^{\prime}={\bf n}$ is the time-dependent coordinate transformation of the unit vector ${\bf n}$ that corresponds to the coordinates $(\alpha,\delta)$ in the right-handed equatorial system. Here, we adopt the convention used by Ahlers et al. (2016) where roman indices ($i$, $j$) refer to pixels in the local sky map and fraktur indices ($\mathfrak{a}$) refer to pixels in the celestial sky map while time bins are indicated by greek indices ($\tau$, $\kappa$). The data observed at a fixed sidereal time bin $\tau$ can be described in terms of the observation in local horizontal sky with bin $i$ as $n_{\tau i}$ or transformed into the celestial sky map with bin $\mathfrak{a}$ as $n_{\tau\mathfrak{a}}$.

The likelihood of observing $n$ cosmic rays is then given by the product of Poisson probabilities

[TABLE]

where $n=\sum_{i,\tau}n_{\tau i}$. We maximize the likelihood ratio of signal over null hypothesis of no anisotropy ($I^{(0)}$, $\mathcal{N}^{(0)}$, $\mathcal{A}^{(0)}$),

[TABLE]

with $I_{\mathfrak{a}}^{(0)}=1$. The maximum likelihood estimators of $\mathcal{A}_{i}$ and $\mathcal{N}_{\tau}$ are then

[TABLE]

given the boundary condition

[TABLE]

In this combined analysis of HAWC and IceCube data, the likelihood (Eq. A3) is generalized to a product over data sets with individual detector exposures but the same relative intensity. The total accumulated exposure $\mathcal{E}$ in Eq. (A1) becomes a sum over disjoint sky sectors, whose union covers the entire field of view. In this analysis the integrated field of view of each detector corresponds to a sector. As before, we assume that the exposure in each sector can be expressed as a product of its angular-integrated exposure $E^{s}$ and relative acceptance in terms of azimuth $\varphi$ and zenith angle $\theta$ as

[TABLE]

The values of $I$, $\mathcal{N}$, and $\mathcal{A}$ of the maximum likelihood ratio (Eq. A4) $(I^{\star},~{}\mathcal{N}^{\star},~{}\mathcal{A}^{\star})$ must obey the implicit relations

[TABLE]

Here, we have introduced the window function $w^{s}_{i}$ of the sector $s$ which is equal to $1$ if the pixel $i$ is located in the sector and [math] otherwise. The binned quantity $\mathcal{A}^{s}_{\tau\mathfrak{a}}$ in Eq. (A9), corresponds to the relative acceptance of sector $s$ seen in the equatorial coordinate system in pixel $\mathfrak{a}$ during time bin $\tau$. Equations (A9) to (A11) correspond to a nonlinear set of equations that cannot be solved in explicit form but one can iteratively approach the best-fit solution.

This reconstruction method is a simple generalization of the iterative method outlined in  Ahlers et al. (2016), where now the relative acceptances $\mathcal{A}$ and isotropic expectation $\mathcal{N}$ for each detector are evaluated as independent quantities. This is a valid approach as long as the rigidity distributions of the data sets are very similar.