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

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 2013; Sveshnikova et al. 2013; Kumar & Eichler 2014; 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, 2017; Guillian et al. 2007; Abdo et al. 2008, 2009; Aglietta et al. 2009; Munakata et al. 2010; Abbasi et al. 2011, 2010, 2012; De Jong 2011; Aartsen et al. 2013, 2016; Bartoli et al. 2013, 2015, 2018; Abeysekara et al. 2014, 2018b). The observed large-scale anisotropy has an amplitude of about 10⁻³ and small-scale structures of amplitude of 10⁻⁴ with angular size of 10°–30°.

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_ℓm bias the interpretation of the cosmic ray distributions in the context of particle diffusion in the local interstellar medium. 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 ℓ (large scale), which is most affected by partial sky coverage. It should be noted that neither observatory is sensitive to variations across decl. bands, as events recorded from a fixed direction in the local coordinate system can only probe the cosmic-ray flux in a fixed decl. band δ. 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.

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°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 m² array of 300 close-packed water Cerenkov 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 PMTs are anchored in an equilateral triangle of side length 3.2 m, 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 Cerenkov 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% uptime 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 and an angular resolution of 04 for energies above 10 TeV, HAWC is an ideal instrument for measuring the cosmic-ray arrival direction distribution with unprecedented precision.

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 5160 optical sensors, called digital optical modules (DOMs), frozen in the ice at depths from 1450 to 2450 m below the surface of the ice. A DOM consists of a pressure-protective glass sphere that houses a 10-inch Hamamatsu PMT 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 Cerenkov 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° according to simulation (Abbasi et al. 2011) but worsens past zenith angles of approximately 70° (Aartsen et al. 2017). This is not to be confused with the ∼06 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 described in Section 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 data set used in most neutrino analyses.

The data set selected for this analysis is composed of 5 years of data collected by the IceCube Neutrino Observatory between 2011 May and 2016 May, as well as 2 years of data from the HAWC Gamma-Ray Observatory collected between 2015 May and 2017 May. 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⁻⁵, which is smaller than the statistical error of this analysis (see Section 7.2). Cuts are applied to each data set 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 θ > 57°, 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_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 }$, which describes the "clumpiness" of the air shower (Abeysekara et al. 2017) with ${ \mathcal P }\gt 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 ζ_i = log₁₀(Q_eff,i). For each PMT hit, i, an expectation is assigned $\langle {\zeta }_{i}\rangle$ by averaging the ζ_i in all PMTs contained in an annulus containing the hit, with a width of 5 m, centered at the core of the air shower. ${ \mathcal P }$ is then calculated using the χ² formula:

$$\begin{eqnarray}&&{ \mathcal P }=\displaystyle \frac{1}{N}\displaystyle \sum _{i=0}^{N}\displaystyle \frac{{\left({\zeta }_{i}-\langle {\zeta }_{i}\rangle \right)}^{2}}{{{\sigma }_{{\zeta }_{i}}}^{2}}.\end{eqnarray} \tag{ 1 }$$

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 (R log L < 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}}\gt 9\cos (\theta )$ and ${l}_{\mathrm{dir}}\gt 200\cos (\theta )$ meters. This cut depends on the reconstructed zenith angle θ 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 Section 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₁₀(E/GeV) (Aartsen et al. 2016).

Table 1.  Comparison of the IceCube and HAWC Data Sets

	IceCube		HAWC
Latitude	90°S		19°N
Detection method	Muons produced by CR		Air showers produced by CR and γ
Field of view	−90°/−16° (δ), ∼4 sr (same sky over 24 hr)		−30°/68° (δ), ∼2 sr (8 sr observed/24 hr)
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°–3°	2°–6°	04–08	04–10
Events	2.8 × 1011	1.7 × 1011	7.1 × 1010	2.8 × 1010

Note. The median primary particle energy, angular resolution, and number of remaining events is shown in the subcolumns after applying only quality cuts and after applying both energy and quality cuts. The angular resolution of IceCube corresponds to the DST data set 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.

Download table as:  ASCIITypeset image

Figure 1 shows the distribution of data as a function of decl. The resulting energy distribution of the two data sets 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–40 TeV, in the case of IceCube, and 2.5–30 TeV for HAWC around the median energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The two experiments have different responses to the cosmic ray mass composition. This is largely due to the detection method. Particles entering Earth's atmosphere (15–20 km above sea level) interact with nuclei in the air and produce a cascade of secondary particles. This particle cascade continues to grow until ionization becomes the dominant energy-loss mechanism. The depth X_max 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, as 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 (Hörandel 2003) and detailed simulations of the detector response.

We compute the relative intensity as a function of J2000 equatorial coordinates (α, δ) by binning the sky into an equal-area grid with a bin size of 09 using the HEALPix library (Górski et al. 2005). The angular distribution can be expressed as ϕ(α, δ) = ϕ^isoI(α, δ), where ϕ^iso corresponds to the isotropic flux (i.e., the flux averaged over the full celestial sphere), and I(α, δ) is the relative intensity of the flux as a function of R.A. α and decl. δ 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 δ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 δ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,

$$\begin{eqnarray}&&\delta {I}_{{\mathfrak{a}}}=\displaystyle \frac{{N}_{{\mathfrak{a}}}-\langle {N}_{{\mathfrak{a}}}\rangle }{\langle {N}_{{\mathfrak{a}}}\rangle }.\end{eqnarray} \tag{ 2 }$$

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⁻³, 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. δ) 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, because 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

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 σ) is then calculated as

$$\begin{eqnarray}&&{S}_{i}=\sqrt{2}{\left(-{\mu }_{i,\mathrm{on}}+{\mu }_{i,\mathrm{off}}+{n}_{i}\mathrm{log}\displaystyle \frac{{\mu }_{i,\mathrm{on}}}{{\mu }_{i,\mathrm{off}}}\right)}^{1/2}.\end{eqnarray} \tag{ 3 }$$

For each pixel i in the celestial sky, we define expected on-source and off-source event counts from neighbor pixels in a disk 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

$$\begin{eqnarray}&&{n}_{i}=\displaystyle \sum _{j\in {{ \mathcal D }}_{i}}\displaystyle \sum _{\tau }{n}_{\tau j},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\mu }_{i,\mathrm{on}}=\displaystyle \sum _{j\in {{ \mathcal D }}_{i}}\displaystyle \sum _{\tau }{{ \mathcal A }}_{\tau j}{{ \mathcal N }}_{\tau }{I}_{j},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\mu }_{i,\mathrm{off}}=\displaystyle \sum _{j\in {{ \mathcal D }}_{i}}\displaystyle \sum _{\tau }{{ \mathcal A }}_{\tau j}{{ \mathcal N }}_{\tau }{I}_{j}^{\mathrm{reference}},\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal A }}_{\tau j}$ is the relative acceptance of the detector in pixel j and sidereal time bin τ, ${{ \mathcal N }}_{\tau }$ gives the expected number of isotropic events in sidereal time bin τ, I_j is the relative intensity, and I = I^reference + I^residual is divided into a contribution from the reference map and the residual relative intensity. For small-scale features, I^reference corresponds to the first three spherical harmonic components (ℓ ≤ 3) of the relative intensity. In order to distinguish excess and deficit, we multiply Equation (3) by the sign of each smoothed pixel δ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^ℓm,

$$\begin{eqnarray}&&\delta I({{\boldsymbol{u}}}_{i})=\displaystyle \sum _{{\ell }=1}^{\infty }\displaystyle \sum _{m=-{\ell }}^{{\ell }}{a}_{{\ell }m}{Y}_{{\ell }m}({{\boldsymbol{u}}}_{i}).\end{eqnarray} \tag{ 7 }$$

The vector components of the dipole in terms of the spherical harmonic expansion Y_ℓm in equatorial coordinates are related to the a_ℓm coefficients with

$$\begin{eqnarray}&&{\boldsymbol{\delta }}\equiv ({\delta }_{0h},{\delta }_{6h},{\delta }_{N})=\sqrt{\displaystyle \frac{3}{2\pi }}(-{\mathfrak{R}}({a}_{11}),{\mathfrak{I}}({a}_{11}),{a}_{10}),\end{eqnarray} \tag{ 8 }$$

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

From Equation (8) and the a_ℓm coefficients, one can obtain the horizontal components of the dipole δ_0h and δ_6h with respect to the 0h and 6h R.A. axes. The phase and amplitude of the projected dipole on the equatorial plane are given by

$$\begin{eqnarray}&&({\delta }_{0h},{\delta }_{6h})=({\tilde{A}}_{1}\cos {\alpha }_{1},{\tilde{A}}_{1}\sin {\alpha }_{1}),\end{eqnarray} \tag{ 9 }$$

where α₁ 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₁ through the dipole inclination δ₀ 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

$$\begin{eqnarray}&&{{ \mathcal C }}_{{\ell }}=\displaystyle \frac{1}{2{\ell }+1}\displaystyle \sum _{m=-{\ell }}^{{\ell }}| {a}_{{\ell }m}{| }^{2},\end{eqnarray} \tag{ 10 }$$

for each value of ℓ. 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 }}$:

$$\begin{eqnarray}&&{\tilde{{ \mathcal C }}}_{{\ell }}=\displaystyle \frac{1}{2{\ell }}\displaystyle \sum _{m=-{\ell },m\ne 0}^{{\ell }}| {a}_{{\ell }m}{| }^{2}.\end{eqnarray} \tag{ 11 }$$

The angular power spectrum provides an estimate of the significance of structures at different angular scales of ∼180°/ℓ. In the ideal case of a 4π sky coverage, the multipole moments a_ℓ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.

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 data set compared to that from HAWC at energies of ∼10 TeV.

Zoom In Zoom Out Reset image size

Figure 5 is the residual small-scale anisotropy after subtracting the fitted multipole from the spherical harmonic expansion with ℓ ≤ 3 from the large-scale map in Figure 4 in order to reveal structures smaller than 60°. 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 data sets. The one referred to as "region A" by the Milagro Collaboration (Abdo et al. 2008) roughly extends from (54°, −16°) to (78°, 18°), in equatorial coordinates (δ, α). 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We obtain the a_ℓ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_ℓm values in Table 3 are δ_0h = 9.16 × 10⁻⁴ and δ_6h = 7.25 × 10⁻⁴ (±0.04 × 10⁻⁴), respectively, with respect to the 0h and 6h R.A. axes. The dipole amplitude and phase ${\tilde{A}}_{1}=(1.17\pm .01)\,\times {10}^{-3}$, ${\alpha }_{1}=38\buildrel{\circ}\over{.} 4\pm 0\buildrel{\circ}\over{.} 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 δα₁ ∼ 26, respectively (see Section 7).

The angular power spectrum for the combined data set in Figure 7 provides an estimate of the significance of structures at different angular scales of ∼180°/ℓ. 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 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 data sets 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), noted previously.

Zoom In Zoom Out Reset image size

The measured quadrupole component has an amplitude of 6.8 × 10⁻⁴ and is inclined at 207 ± 03 above (and below) the 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₂₁ and a₂₂ nonvertical 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 (∼10⁻³). 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 ℓ = 2 and ℓ = 3 multipole components and reduces correlations between spherical harmonic expansion coefficients a_ℓm.

7.1. Overlapping Region

We have studied two adjacent δ 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 data sets 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 data sets (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 α = 50° 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 the increasing zenith angle, as is the quality of angular reconstructions. As a result, δ bins closer to the horizon contain a high level of contamination from bins in higher zenith angles.

Zoom In Zoom Out Reset image size

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₁₁ and a₁₋₁, 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 ℓ pseudo-modes under partial sky coverage that primarily affects the multipolar components ℓ = 2, ℓ = 3, and to a lesser degree, ℓ = 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 δ_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. δ_max, keeping δ_min = −90°.

Zoom In Zoom Out Reset image size

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⁻⁵ in this analysis. Figure 10 shows the correlation matrix (Efstathiou 2004) of the different ℓ-modes up to ℓ = 30, calculated using the PolSpice⁷⁵ software package. The correlation between ℓ-modes due to partial sky coverage is appreciable for larger ℓ, though to a lesser degree.

Zoom In Zoom Out Reset image size

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, which 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 represent a systematic effect of the solar frame on the sidereal anisotropy (Aartsen et al. 2016). The anti-sidereal distribution of the HAWC data set has a maximum amplitude of 5 × 10⁻⁵. 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 ∼3 × 10⁻⁵. The worst-case uncertainty on the reconstructed phase of the dipole is δα = 26 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

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}I}{I}=(\gamma +2)\displaystyle \frac{v}{c}\cos ({\theta }_{v}),\end{eqnarray} \tag{ 12 }$$

where I is the cosmic-ray intensity, γ 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 θ_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° such that after 1 year, the effect is ideally completely canceled 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⁻⁴ 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 ∼10⁻⁵ for the HAWC data set, 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 X_max 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₀ = p/g, where g = 9.87 m s⁻² 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 ${\rm{\Delta }}\{\mathrm{ln}R\}$ and the surface pressure change ΔP is

$$\begin{eqnarray}&&{\rm{\Delta }}\{\mathrm{ln}R\}=\beta \cdot {\rm{\Delta }}P,\end{eqnarray} \tag{ 13 }$$

where β 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 Equation (13) in order to weight individual events. This yields a barometric coefficient of β = −0.0071 hPa⁻¹. The residual contamination from atmospheric variations is estimated to be on the order of ∼10⁻⁶. Temperature variations in the stratosphere can introduce a similar effect with a 24-hour 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, as 1 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 × 10⁻⁴ 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.

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 information regarding 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

8.2. Large-scale Anisotropy and the LIMF

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 (∼115° R.A.) is shown in Figure 12. The fit yields a vector pointed toward (α_fit, δ_fit) = (2292 ± 35, 114 ± 30) in J2000 equatorial coordinates, as shown in Figure 11, along with the corresponding equator (the crossed black curve). The direction is located 9° 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 65 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 α obtained from the dipole fit and the value of δ obtained from the quadrupole fit. The errors on the fit are derived from the χ² distribution shown in Figure 12 and do not include possible systematics uncertainties from the missing m = 0 dipole component.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Magnetic Field Alignment

Source	R.A. [°]	Decl. [°]	Δψ [°]	δN [10−4]
Funsten et al. (2013)	229.7 ± 1.9	23.3 ± 2.5	9.0	−5.03
Frisch et al. (2015)	237.9 ± 16	22.2 ± 16	12.2	−4.77
Zirnstein et al. (2016)	234.4 ± 0.7	16.3 ± 0.6	6.5	−3.42
Boundary Fit	229.2 ± 3.5	11.4 ± 3.0	⋯	−2.36
Dipole/Quadrupole	218.4 ± 0.3 (±2.6)	20.7 ± 0.3 (±2.6)	⋯	−4.41

Note. 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 quadrupole component. The second to last column corresponds to the angular distance Δψ 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 toward the given decl. Error in parentheses for dipole and quadruple correspond to systematic uncertainties.

Download table as:  ASCIITypeset image

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, as we expect diffusion to be anisotropic with the fastest propagation along the magnetic field lines (Effenberger et al. 2012; Kumar & Eichler 2014; 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₁ through the dipole inclination δ₀ 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}_{-2.0}^{+1.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.

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° to +76° 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 (α_fit, δ_fit) = (2292 ± 35, 114 ± 30) 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}$, α₁ = 384 ± 03. Based on the assumption that the true dipole is aligned along the LIMF, we estimated the missing vertical component to be δ_N  ∼ $-{3.97}_{-2.0}^{+1.0}\times {10}^{-4}$.

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 support from USA—US 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), US 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).

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, Robitaille et al. 2013; Price-Whelan 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).

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{\rm{\Omega }}A({\rm{\Omega }})=1$),

$$\begin{eqnarray}&&{ \mathcal E }(t,\varphi ,\theta )\simeq E(t)A(\varphi ,\theta ),\end{eqnarray} \tag{ 14 }$$

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 ΔΩ_i of the local coordinate sphere corresponding to coordinates (θ_i, φ_i) in a sidereal time interval Δt_τ is

$$\begin{eqnarray}&&{\mu }_{\tau i}\simeq {I}_{\tau i}{{ \mathcal N }}_{\tau }{{ \mathcal A }}_{i},\end{eqnarray} \tag{ 15 }$$

where ${{ \mathcal N }}_{\tau }\equiv {\rm{\Delta }}{t}_{\tau }{\phi }^{\mathrm{iso}}E({t}_{\tau })$ gives the expected number of isotropic events in sidereal time bin τ independent of pixel, ${{ \mathcal A }}_{i}$ is the relative acceptance of the detector for pixel i, and ${I}_{\tau i}\equiv I({\boldsymbol{R}}({t}_{\tau }){\boldsymbol{n}}^{\prime} ({{\rm{\Omega }}}_{{\mathfrak{i}}}))$ is the relative intensity observed in local coordinates during time bin τ. R(t) n' = n is the time-dependent coordinate transformation of the unit vector n that corresponds to the coordinates (α, δ) 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 (τ, κ). The data observed at a fixed sidereal time bin τ can be described in terms of the observation in the local horizontal sky with bin i as n_τ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,

$$\begin{eqnarray}&&{ \mathcal L }(n| I,{ \mathcal N },{ \mathcal A })=\displaystyle \prod _{\tau i}\displaystyle \frac{{\left({\mu }_{\tau i}\right)}^{{n}_{\tau i}}{e}^{-{\mu }_{\tau i}}}{{n}_{\tau i}!},\end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{ \mathcal L }(n| I,{ \mathcal N },{ \mathcal A })}{{ \mathcal L }(n| {I}^{(0)},{{ \mathcal N }}^{(0)},{{ \mathcal A }}^{(0)})},\end{eqnarray} \tag{ 17 }$$

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

$$\begin{eqnarray}&&{{ \mathcal N }}_{\tau }^{(0)}=\displaystyle \sum _{i}{n}_{\tau i},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{ \mathcal A }}_{i}^{(0)}=\displaystyle \sum _{\tau }{n}_{\tau i}\Space{0ex}{3.0ex}{0ex}/\displaystyle \sum _{\kappa j}{n}_{\kappa j},\end{eqnarray} \tag{ 19 }$$

given the boundary condition

$$\begin{eqnarray}&&\displaystyle \sum _{i}{{ \mathcal A }}_{i}=1.\end{eqnarray} \tag{ 20 }$$

In this combined analysis of HAWC and IceCube data, the likelihood (Equation (16)) is generalized to a product over data sets with individual detector exposures but the same relative intensity. The total accumulated exposure ${ \mathcal E }$ in Equation (14) 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 φ and zenith angle θ as

$$\begin{eqnarray}&&{ \mathcal E }(t,\varphi ,\theta )\simeq \displaystyle \sum _{\mathrm{sector}\,s}{E}^{s}(t){{ \mathcal A }}^{s}(\varphi ,\theta ).\end{eqnarray} \tag{ 21 }$$

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

$$\begin{eqnarray}&&{I}_{{\mathfrak{a}}}^{\star }=\displaystyle \sum _{\tau }{n}_{\tau {\mathfrak{a}}}\Space{0ex}{2.5ex}{0ex}\displaystyle \sum _{s\kappa }{{ \mathcal A }}_{\kappa {\mathfrak{a}}}^{s\,\star }{{ \mathcal N }}_{\kappa }^{s\,\star },\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{ \mathcal N }}_{\tau }^{s\,\star }=\displaystyle \sum _{i}{w}_{i}^{s}{n}_{\tau i}\Space{0ex}{2.5ex}{0ex}/\displaystyle \sum _{j}{{ \mathcal A }}_{j}^{s\,\star }{I}_{\tau j}^{\star },\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{{ \mathcal A }}_{i}^{s\,\star }=\displaystyle \sum _{\tau }{w}_{i}^{s}{n}_{\tau i}\Space{0ex}{3.0ex}{0ex}/\displaystyle \sum _{\kappa }{{ \mathcal N }}_{\kappa }^{s\,\star }{I}_{\kappa i}^{\star }.\end{eqnarray} \tag{ 24 }$$

Here, we have introduced the window function ${w}_{i}^{s}$ of the sector s, which is equal to 1 if the pixel i is located in the sector and 0 otherwise. The binned quantity ${{ \mathcal A }}_{\tau {\mathfrak{a}}}^{s}$ in Equation (22) corresponds to the relative acceptance of sector s seen in the equatorial coordinate system in pixel ${\mathfrak{a}}$ during time bin τ. Equations (22)–(24) 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.