{'fraction_non_alphanumeric': 0.03399525145693935, 'fraction_numerical': 0.017294409669760414, 'mean_word_length': 4.535394265232975, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A New Maximum-Likelihood Technique for Reconstructing Cosmic-Ray Anisotropy at All Angular Scales

The arrival directions of TeV-PeV cosmic rays show weak but significant anisotropies with relative intensities at the level of one per mille. Due to the smallness of the anisotropies, quantitative studies require careful disentanglement of detector effects from the observation. We discuss an iterative maximum-likelihood reconstruction that simultaneously fits cosmic ray anisotropies and detector acceptance. The method does not rely on detector simulations and provides an optimal anisotropy reconstruction for ground-based cosmic ray observatories located in the middle latitudes. It is particularly well suited to the recovery of the dipole anisotropy, which is a crucial observable for the study of cosmic ray diffusion in our Galaxy. We also provide general analysis methods for recovering large- and small-scale anisotropies that take into account systematic effects of the observation by ground-based detectors.


INTRODUCTION
During the past decade, a number of cosmic-ray, γray, and neutrino observatories have found anisotropies in the arrival directions of Galactic cosmic rays at TeV and PeV energies (Guillian et al. 2007;Amenomori et al. 2006;Abdo et al. 2008Abdo et al. , 2009Abbasi et al. 2011Abbasi et al. , 2012Amenomori et al. 2012;Aartsen et al. 2014;Bartoli et al. 2013;Abeysekara et al. 2014; Bartoli et al. 2015). For a recent summary of experimental results, we refer to Di Sciascio & Iuppa (2014). The statistics and resolution of these experiments allow for a two-dimensional reconstruction in the form of anisotropy sky maps. The studies have revealed significant anisotropies at both large and small angular scales. At the largest scales, the anisotropy is approximately dipolar and has a relative intensity on the order of 10 −3 . An explanation of the amplitude and phase of the dipole anisotropy is challenging, but the observations are qualitatively consistent with diffusive propagation of cosmic rays from Galactic sources (Ptuskin et al. 2006;Erlykin & Wolfendale 2006;Blasi & Amato 2012;Mertsch & Funk 2015).
Medium-and small-scale anisotropies have been observed at the 10 −4 level. These features are less understood but could be a combined effect of nearby cosmicray sources (Salvati & Sacco 2008;Biermann et al. 2013) and the local interstellar magnetic field structure, which can introduce an energy-dependent magnetic mirror leakage (Drury & Aharonian 2008), preferred cosmic-ray transport directions (Malkov et al. 2010), or magnetic lenses (Battaner et al. 2011(Battaner et al. , 2015. The observed power spectrum of cosmic-ray anisotropies agrees well with the expected effect of cosmic-ray scattering in the local magnetic field (Giacinti & Sigl 2012;Ahlers 2014;Ahlers & Mertsch 2015;López-Barquero et al. 2015). The possible influence of the heliosphere via magnetic reconnections in the heliotail (Lazarian & Desiati 2010), nonisotropic particle transport in the heliosheath (Desiati & Lazarian 2013) or the heliospheric electric field structure (Drury 2013) has also been considered. More exotic origin models invoke strangelet production in molecular clouds (Kotera et al. 2013) or in neutron stars (Perez-Garcia et al. 2014).
The measurement of cosmic-ray anisotropies with relative intensity below 10 −3 is an observational challenge, since it is necessary to account for minuscule variations in the acceptance and uptime of the detector carrying out the measurement. For illustration, Fig. 1 shows a simulated realization of the cosmic-ray anisotropy following the model of Ahlers (2014). This model was chosen due to its qualitative agreement with the observation of anisotropy in TeV cosmic rays by IceCube (Aartsen et al. 2014) and HAWC (Abeysekara et al. 2014). The cosmic-ray anisotropy is expected to remain constant in the celestial coordinate system over the period of the observation. Ground-based cosmic observatories with a limited field of view, however, are exposed to different parts of the celestial sphere as the Earth rotates during one sidereal day. As an example, the sky map of Fig. 1 indicates the instantaneous field of view of the HAWC detector (Abeysekara et al. 2014) (at latitude 19 • N) at a local sidereal time of 9h. Hence, the observed event distribution accumulated over many sidereal days depends not only on the cosmic-ray anisotropy but also on the nonuniform and time-dependent detector exposure. It is therefore necessary to reconstruct a reference map that represents the response of the detector to an isotropic cosmic-ray flux. For ground-based detectors, in which the atmosphere acts as part of the observatory, sufficiently accurate simulations of the detector exposure are usually not achievable. Therefore, most analyses use the data themselves to estimate the relative intensity and detector exposure simultaneously.
Examples of this technique are the timescrambling (Alexandreas et al. 1993) and directintegration (Atkins et al. 2003)  rate of events observed in a detector as a function of local sidereal time is integrated against the relative acceptance of the detector during an integration period (or scrambling interval) ∆t. The idea behind this method is that variations in the event rates introduced by the cosmic-ray anisotropies will average out to some extent as the detector observes different parts of the celestial sphere over the course of each sidereal day. The result is an estimate of the number of events expected in the detector between time t and t + ∆t. Subtracting the expected number of events from the actual observations yields a residual counts map which can be explored for anisotropy. The integration time interval ∆t acts as an effective smoothing parameter for the counts map, since the method averages cosmic-ray arrival directions over angular scales of 15 • (∆t/1hr). In principle, choosing ∆t = 24h would produce a residual counts map with features covering the full sky (360 • ). However, for detectors located in the middle latitudes, the instantaneous exposure of the detector does not match the full daily exposure, since the 24h integrated field of view is much larger than the instantaneous field of view, cf. Fig 1. As a result, large-scale structures in the residual counts map, in particular the dipole, are strongly attenuated when using these methods.
To improve estimates of large-scale anisotropy using detectors in the middle latitudes, we describe a maximum likelihood construction that can be used to disentangle the anisotropy from detector effects. The technique is based on the same ansatz used by the time-scrambling or direct-integration method, that the total accumulated exposure of the detector can be factorized into a timedependent event rate and a time-independent relative acceptance map. We begin by describing the technique in Section 2. In Section 3, we apply the maximumlikelihood method to simulated data and show that largeand small-scale anisotropies can be reconstructed with relatively little of the distortion observed in the directintegration or time-scrambling techniques. We compare our method to alternative techniques in Section 4. We then discuss analysis methods of large-and small-scale anisotropies of the reconstructed anisotropy maps in Sec-tion 5 before concluding in Section 6.

MAXIMUM LIKELIHOOD METHOD
In the following, we will assume that the total accumulated detector exposure E can be expressed as a product of its angular-integrated exposure E and relative acceptance A in terms of azimuth angle ϕ (from north increasing to the east) and zenith angle θ as (1) Without loss of generality, we require that the relative acceptance is normalized to dΩA(ϕ, θ) = 1. This approximation assumes that the relative acceptance of the detector remains approximately constant over time. The ansatz is identical to the approach used in direct integration or time scrambling. Let us also assume that the flux of cosmic rays at the energies of interest remains constant as a function of time, varying only as a function of celestial longitude α (right ascension) and latitude δ (declination). Due to the strong diffusion of cosmic rays in the Galactic environment, the flux is dominated by an isotropic term φ iso . Hence, the total flux can be expressed as where I(α, δ) is the relative intensity of the flux as a function of position in the sky. The anisotropy is defined as the deviation δI = I − 1 1. Note that this ansatz ignores anisotropies associated with the relative motion of the Earth with respect to the Sun. We will come back to this subtlety in the discussion section.
The local horizontal coordinate system and the celestial (or equatorial) coordinate system are related via a time-dependent transformation. We define n = (cos α cos δ, sin α cos δ, sin δ) as the unit vector corresponding to the coordinates (α, δ) in the right-handed equatorial system. Similarly, the unit vector corresponding to the coordinates (θ, ϕ) in the right-handed local system is n = (cos ϕ sin θ, − sin ϕ sin θ, cos θ). The two unit vectors are related via a time-dependent coordinate transformation n = R(t)n. For an experiment located at geographic latitude Φ and longitude Λ (measured east from Greenwich), the transformation is where ω = 2π/24h and the local sidereal time t is related to the sidereal time at Greenwich t by t = t + Λ/ω.
To simplify calculations on the local and celestial spheres, the sky is binned into pixels of equal area ∆Ω using the HEALPix parametrization of the unit sphere (Gorski et al. 2005). To make the equations more transparent, we use roman indices for pixels in the local sky map and fraktur indices for pixels in the celestial sky map. Time bins are indicated by greek indices. For instance, the data observed at a fixed sidereal time bin τ can be described in terms of the observation in local horizontal sky with bin i as n τ i or transformed into the celestial sky map with bin a as n τ a .
Consider an angular element of the local coordinate sphere ∆Ω i corresponding to coordinates (θ i , ϕ i ). The number of cosmic rays expected from this location in a sidereal time interval ∆t τ with central value t τ is where N τ ≡ ∆t τ φ iso E(t τ ) gives the expected number of isotropic background events in sidereal time bin τ . The quantity A i ≡ ∆Ω i A(θ i , ϕ i ) is the binned relative acceptance of the detector for angular element i, and I τ i ≡ I(R(t τ )n (Ω i )) is the relative intensity observed in the local horizontal system during time bin τ . For simplicity, all expressions that follow assume equal bin sizes ∆Ω = 4π/N pix on the local and celestial spheres and N time equal sidereal bins with ∆t = 24h/N time . However, a more general binning is also possible. Given µ τ i , the likelihood of observing n cosmic rays is given by the product of Poisson probabilities where n τ i is the number of events observed in the local pixel i during time bin τ . This likelihood can be maximized to provide estimators of the relative acceptance function A i and the expected isotropic background count N τ . Consider the case of the null hypothesis of no anisotropy, which we denote I (0) a = 1. Given the boundary condition i A i = 1, the maximum likelihood estimators of A i and N τ are To allow for the possibility of anisotropy, we maximize the likelihood ratio of signal over null hypothesis in N , A, and I.

Invariance under Declination-Scaling
We will demonstrate in the following that the maximum-likelihood method can be used to search for anisotropy on all angular scales. However, before we continue, note that events recorded in a fixed position in the local coordinate system can only probe the cosmic-ray flux in a fixed declination band δ. Hence, the expectation values (4) are invariant under the rescaling where a(δ) is an arbitrary function of declination and the normalization factors b and c are defined such that τ A τ = 1 and a δI a = 0 for the new values. In other words, the maximum-likelihood method is sensitive to anisotropy in right ascension but is insensitive to variations in intensity across declination bands. This is also known to be a limitation of other reconstruction methods, like direct integration or time scrambling ( Because of this degeneracy, we must choose whether to account for the local excess of cosmic rays as originating in an anisotropy in their relative intensity or originating in a variation in the local acceptance. A natural choice is that the anisotropy is normalized to dαδI(α, δ) = 0 for all declinations δ, consistent with the definition dΩδI(α, δ) = 0. This condition can also be formulated in terms of a spherical harmonics expansion of the relative intensity in the equatorial coordinate system (α, δ), as pointed out by Iuppa & Di Sciascio (2013). In general, the relative intensity can be decomposed as a sum over spherical harmonics Y m , Our normalization condition can then be expressed as the condition a 0 = 0 for all . This projection significantly reduces the reconstruction of the low-multipole components of the anisotropy, as we will discuss in the following sections. Note that the true multipole moments a m are an (infinite) superposition of the pseudo multipole moments a m , which are defined as in Eq. (12), but for the product of the relative intensity with the weight function w of the field of view. Provided that the weight function is azimuthally symmetric, w(α, δ) = w(δ), the true multipole moments a 0 are a linear superposition of pseudo multipole moments a 0 . In practice, we can hence use the normalization condition a 0 = 0 for all to ensure a 0 = 0 for all . In terms of the binned relative intensity and weight function this is equivalent to the condition a w a Y 0 a δI a = 0 for all .

Maximum Likelihood Algorithm
The maximum (I , N , A ) of the likelihood ratio (8) must obey the implicit relations together with a w a Y 0 a δI a = 0 and i A i = 1. In Eq. (13), we introduced the binned quantity A τ a ≡ ∆Ω a A(R T (t τ )n(Ω a )) corresponding to the relative acceptance seen in the equatorial coordinate system in pixel a during time bin τ .
Equations (13) to (15) correspond to a nonlinear set of equations that cannot be solved in an explicit form. However, one can approach the best-fit solution via the following iterative method: (i) Initialize at the maximum of the null hypothesis, (I (0) , N (0) , A (0) ).
Note that the cosmic-ray anisotropy obtained in the first iteration step, δI (1) , corresponds to the result that would be obtained by the method of direct integration (Atkins et al. 2003). This method was also used in the recent HAWC analysis (Abeysekara et al. 2014). The successive iteration steps of the maximum-likelihood method reoptimize the relative acceptance A and the isotropic background rate N for a given anisotropy. We will study this optimization process more quantitatively in the following section.

SIMULATION AND PERFORMANCE
To demonstrate the reconstruction of a cosmic-ray anisotropy with the maximum likelihood technique, we simulated a set of arrival directions based on the anisotropy shown in Fig. 1. The simulation is based on a random realization of a relative intensity δI, which follows a power-law spectrum with ≥ 1 of the form C = 10 −7 (18/(2 + 1)/( + 1)/( + 2)), after the model of Ahlers (2014). This model was chosen due to its qualitative agreement with the observation of anisotropy in TeV cosmic rays by IceCube (Aartsen et al. 2014) and HAWC (Abeysekara et al. 2014). The celestial sphere is binned following the HEALPix parametrization (Gorski et al. 2005) with parameter n side = 64 corresponding to N pix = 49, 152 pixels with a binsize of about 1 • diameter.
The simulated detector is located at geographic coordinates (Φ = 19 • N and Λ = 97 • W), the location of the HAWC observatory (Abeysekara et al. 2014). The instantaneous field of view is restricted to zenith angles below 60 • . A projection of the instantaneous field of view at 9h local sidereal time onto the equatorial coordinate system was already shown in Fig. 1. We assume that the relative detector acceptance at any time follows A(θ, ϕ) ∝ cos θ[1 + A sin θ sin 2 (ϕ − ϕ 0 )] with A = 0.2 and ϕ 0 = 10 • . The local acceptance maps are reconstructed with the same resolution as the anisotropy map.
The expected isotropic event number is binned in N time = 360 sidereal time bins with 4 min bin size and normalized such that τ N τ 5 × 10 10 , corresponding to the integrated event number of the recent HAWC analysis (Abeysekara et al. 2014). We introduce a simple statistical toy model for the variation of the expected background rate N τ ; the accumulated data are expected to be stable over periods sampled from an exponential distribution with expectation value of 20 min (5 sidereal time bins). In each stable period, we assume that N τ has a normal fluctuation of 1% around the background expectation. This model mimics the data variation of the analysis by Abeysekara et al. (2014).
The results of the reconstruction 6 are shown in Figs. 2 and 3. The top panel of Fig. 2 shows the optimization of the reconstructed background expectation N (n) as the relative quantity N (n) /N (0) − 1. The iteration renormalizes N (n) at the level of 10 −4 depending on the local sidereal time. The relative features introduced by the optimization process are easy to understand. For instance, the bump at 9h local sidereal time corresponds to the instantaneous field of view indicated in Fig. 1. In comparison with the simulated anisotropy shown in Fig. 1, one can notice that at this time the detector observes a part of the sky with a strong underfluctuation. The initial estimate N (0) is therefore too low and the iteration compensates for this effect. The bottom panel of Fig. 2 shows the consecutive log-likelihood values of the iteration. Note that for this simulated data set the iteration only requires about 10 steps to converge. We run the reconstruction for 10 more iteration steps to verify the convergence and stability of the method.
The map in Fig. 3a shows the expected relative intensity smoothed with a top-hat function with 10 • radius after application of the m = 0 filter. The maps  in Figs. 3b and 3d show the reconstructed anisotropy in the first iteration step (corresponding to the result from direct integration) and after 20 iterations, respectively. Comparing the maps in Figs. 3b and 3a one notices that, qualitatively, the expected small-scale features are already reproduced in the first iteration step. However, the difference map in Fig. 3c indicates, that the residual map has large-scale features that are misreconstructed. On the other hand, the corresponding difference map in Fig. 3e after 20 iterations is closer to the expected map in Fig. 3a.
The residual anisotropy map shown in Fig. 3e is re-lated to the Poisson variation of the event rate. We can make this statement more quantitative via a power spectrum analysis of expected and reconstructed anisotropy maps. The relative intensity can be decomposed as a sum over spherical harmonics, as in Eq. (12). Unfortunately, due to the limited integrated field of view, the true coefficients a m cannot be unambiguously reconstructed. However, for the present discussion of residual anisotropies in the iterative method, it is sufficient to study the pseudo multipole moments, a m , corresponding to the harmonic expansion of the anisotropy multiplied by the weight function of the field of view. In our  Note that the small bias of the reconstructed pseudo power toward higher values for ≥ 6 is due to the noise level of the maps as we will show in Section 5. example, the weight function is simply equal to 1 for declinations −41 • < δ < 79 • and 0 otherwise. From these pseudo multipole moments, we can compute the pseudo power spectrum, Figure 4 shows the pseudo power spectrum of the anisotropy maps for the first and last iteration step in comparison to the expected spectrum corresponding to the true anisotropy with a 0 = 0 and multiplied by the weight function. To estimate the variance of the reconstructed power spectra due to Poisson statistics, we repeat the analysis 100 times with low-resolution maps (n side = 16 with bin size ∆θ 4 • ) and show the central 68% range of the data. One can see that the first iteration step drastically underestimates the power of the dipole ( = 1) and quadrupole ( = 2). This was already qualitatively visible in the map of Fig. 3c. On the other hand, the pseudo power spectrum of the last iteration step agrees well with the low-multipoles. One can also notice a small bias of the ≥ 6 multipoles towards larger values. This is due to the noise level of the maps as we will discuss in detail in Section 5.

COMPARISON TO OTHER METHODS
Various reconstruction methods have been developed previously to extract the anisotropy from the isotropic background. We have already mentioned the methods of direct integration Atkins et al. (2003) and time scrambling (Alexandreas et al. 1993) in the introduction, which have been used in cosmic-ray anisotropy studies of Super-Kamiokande (Guillian et al. 2007), Tibet-ASγ (Amenomori et al. 2006 The time-integration method is closely related to the maximum-likelihood method presented here. The opti-mal relative detector acceptance (7) and the isotropic background expectation (6) of the null hypothesis (δI = 0) are equivalent to the estimates via direct integration with an integration period of ∆t = 24 h. The map of Fig. 3b shows the resulting anisotropy estimate of this method.
Another two-dimensional reconstruction technique is the Forward-Backward method used by Milagro (Abdo et al. 2009), which is closely related to the onedimensional East-West method (Bonino et al. 2011). These methods analyze the relative right-ascension derivative of event rates, ∂ α n/n, at each sidereal time, either for individual declination bands (in the case of Forward-Backward) or the entire field of view (in the case of East-West). The anisotropy can then be reconstructed from the first derivative, noting that ∂ α δI ∂ α n/n, up to an overall normalization constant in each declination band. This again reflects the invariance (9) to (11) and the inability to reconstruct the a 0 moments in equatorial coordinates.
In general, the Forward-Backward (and the East-West) method has the advantage that the quantity ∂ α n/n does not explicitly depend on the instantaneous background level. For instance, the local deficit of background events at a local sidereal time of 9h in the example shown in Fig. 1 does not affect the estimate of the derivative ∂ α n/n. However, the reconstruction methods discussed in Abdo et al. (2009) assume that the relative detector acceptance is quasi-symmetric under Forward-Backward (or East-West) reflection, corresponding to the transformation (θ, ϕ) → (θ, −ϕ). This condition can be expressed in terms of the asymmetry ≡ [A(θ, ϕ) − A(θ, −ϕ)]/[A(θ, ϕ) + A(θ, −ϕ)] as | | 1. After proper renormalization of the reconstructed derivative terms (Abdo et al. 2009) the leading order effect is an O( 2 ) correction of the anisotropy amplitude. We expect that the condition | | 1 is met by most observatories. However, our method is also applicable for large asymmetries approaching 2 1.
Iterative methods for the reconstruction of the anisotropy have also been developed for and applied to data of Tibet-ASγ (Amenomori et al. 2005(Amenomori et al. , 2010(Amenomori et al. , 2012 and ARGO-YBJ (Bartoli et al. 2015). The Equi-Zenith Angle method (Cui & Yan 2003) estimates the isotropic cosmic-ray background at a celestial bin a to be the average of those acceptance-corrected events that arrived from the same zenith angle band as the field of view wraps around the celestial equator. This corresponds to a generalization of the ansatz (1) to E = s E s A s , where the sum runs over the different equi-zenith-angle sectors with individual background rates E s and relative acceptance A s . It is straightforward to generalize our likelihood-based method with this ansatz (see Appendix A). The likelihood-based iteration method presented in this paper has the advantage that it can be derived from a firm statistical approach. If the relative detector acceptance can be regarded as stable in local sidereal time over the entire field of view, it provides an even simpler iterative reconstruction method.

HARMONIC ANALYSIS
We now turn to the harmonic analysis of the anisotropy. We have already introduced the harmonic expansion of the relative intensity in Eq. (12). Of par-ticular relevance for the theory of cosmic-ray diffusion is the strength of the dipole components ( = 1), but note that TeV-PeV cosmic-ray data also show significant multipole moments at smaller angular scales, like quadrupole ( = 2), octupole ( = 3), etc.
We have already shown in Section 3 that traditional anisotropy methods can significantly underestimate the low-pseudo power spectrum of the data. This effect can be compensated by the iterative method presented in Section 2. However, this method will not compensate for the loss in power of the m = 0 coefficients of the analysis. The observed spectrum is hence always a systematic underestimation of the true anisotropy.
An additional uncertainty comes from the limited integrated field of view of most observatories. 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 a 0 ). However, as already discussed earlier, the partial sky coverage of individual experiments does only allow reconstructing the pseudo multipole moments spectrum of the (reduced) anisotropy multiplied by the weight function.
The pseudo multipole moments a m are related to the true multipole moments a m via a linear transformation (see, e.g., the review by Efstathiou (2004)) where the coupling matrix K depends on the multipole spectrum b m of the weight function of the field of view (see Appendix B). For an azimuthally symmetric weight function the strength of the mixing between moments a m and a m is determined by the moments b k0 with | − | ≤ k ≤ + . This mixing can be small for individual moments as pointed out by Denton & Weiler (2015). However, it is important to emphasize that the full multipole spectrum cannot be unambiguously reconstructed from a partial sky coverage since the infinite-dimensional matrix K is not invertible.
5.1. Large-Scale Anisotropy Whereas the full transition matrix K cannot be inverted, we can attempt an approximate reconstruction of the low-anisotropy via a truncation of the multipole expansion after a maximum . The corresponding truncated matrix K is then invertible. For instance, assuming a pure dipole anisotropy, ≤ 1, and a uniform sky coverage between declination δ 1 and δ 2 gives the dipole transition elements and K 1-11-1 = K 1111 . The dipole strength perpendicular to the Earth's rotation axis can then be estimated as For instance, in our MC simulation with a zenith cut of 60 • and a latitude Φ 19 • , we have δ 1 = −41 • and δ 2 =  79 • giving K 1111 0.92. Hence, the pseudo multipole moments of the projected dipole, a 11 and a 1-1 , have to be corrected by a moderate factor 1/K 1111 1.09 to recover the true momentsâ 11 andâ 1-1 . On the other hand, using the same zenith cut at the location of IceCube/IceTop gives δ 1 = −90 • and δ 2 = −30 • , yielding K 1111 0.16 and a correction factor 1/K 1111 6.4.
However, we emphasize that this treatment is only correct under the assumption that the true anisotropy is dominated by a dipole. In Fig. 5a, we show a map of the reconstructed large-scale anisotropy including dipole ( = 1), quadrupole ( = 2), and octupole ( = 3) with a truncation of the matrix in Eq. (17) after = 3. Note that the model of Ahlers (2014) used for this simulation also assumes significant power at higher multipoles with > 3. The truncation after = 3 is hence not per se justified by this model. However, we can still use the reconstructed large-scale anisotropy map in Fig. 5a as a background model to define the small-scale anisotropy in the full anisotropy map via subtraction.

Small-Scale Anisotropy
The statistical significance of residual anisotropy features in the final reconstructed map is usually estimated using the method introduced by Li & Ma (1983) for applications in gamma-ray astronomy. A direct application of this method does not account for the optimization process of the time-dependent exposure. However, it is rather straightforward to generalize the method of Li & Ma (1983) to our case.
We begin by dividing the reconstructed relative intensity into a contribution of large-scale features and smallscale features, I = I large + I small . For each pixel a in the celestial sky, we define expected on-source and off-source event counts in a disc of radius ψ centered on that pixel. Given the set of pixels D a , the observed and expected counts are The significance map (in units of Gaussian σ) is then calculated as S a = √ 2 −µ a,on + µ a,off + n a log µ a,on µ a,off In Fig. 5b, we show the residual anisotropy map after subtraction of the large-scale anisotropy map shown in Fig. 5a. The corresponding significance map using (24) is shown in Fig. 5c. Usually, the significance is multiplied by the sign of the (smoothed) anisotropy (middle panel) to distinguish deficits (blue) from excesses (red), and we follow here the same convention. We can see that after the subtraction of large-scale features there is still significant power in small-scale features at a significance level of 8σ. Another approach to study small-scale structure is via the power spectrum (16), which quantifies the absolute amplitude of the multipole components but ignores their phases. We have already used the power spectrum to quantify the performance of our iterative method described in Section 2. We can also use this quantity in situations where the small-scale anisotropy structure is dominated by random scattering of cosmic rays in local magnetic fields (Giacinti & Sigl 2012;Ahlers 2014;Ahlers & Mertsch 2015;López-Barquero et al. 2015).
As mentioned in the previous section, the true multipole spectrum cannot be recovered unambiguously from a partial sky coverage, and the same is true for the power spectrum. However, in certain situations one can make additional assumptions about the ensemble-averaged expectation values of the multipole components. In the following, we will assume that the harmonic coefficients are Gaussian random fields, which in the ensemble-average follow a m a * op = δ mp δ o C . In this particular case, we can recover the ensemble-averaged power spectrum C via the relation The transfer matrix M is known from the study of temperature anisotropies in the cosmic microwave background (Efstathiou 2004). However, for our situation of cosmic-ray anisotropies we have again to account for the fact that the m = 0 moments are filtered out by the reconstruction. This leads to a modified expression for M that we discuss and provide in Appendix B. In Eq. (25), we have also introduced the noise power spectrum N , which can in general be calculated from the relative intensity variance from the likelihood function (5). In the following, we will use an approximation that only depends on pixel-by-pixel Poisson noise, which gives a flat spectrum N = N with In Fig. 6, we show the estimated power spectrum inferred from the relative intensity map of the last iteration step shown in Fig. 3d. The dashed line corre-sponds to the expected input power from which the input power spectrum (black data points) is sampled as a Gaussian random field. The sampling introduces a scatter of (∆C ) 2 = 2C 2 /(2 + 1) around the expected input power (cosmic variance). The green data points show the pseudo power spectrum, not accounting for the weight function. The horizontal blue dotted line shows the noise level of Eq. (26). Note that the noise level scales as N ∝ 1/N tot with the total number of events. The red data points are the best estimators of the true power spectrum C . The variance of the pseudo and estimated power spectra are given in Appendix B.
6. SUMMARY In this paper, we have discussed a novel twodimensional cosmic-ray reconstruction method. It is based on a maximum likelihood analysis that provides implicit best-fit expressions for the relative intensity, relative detector acceptance, and background expectations. We have provided a detailed iterative method onhow the relative intensity can be reconstructed from these implicit maximum likelihood solutions. The performance of this likelihood-based method was studied via a simulated example, mimicking the position and performance of the HAWC observatory.
In general, ground-based observatories are insensitive to cosmic-ray anisotropy variations across declination bands. In terms of a spherical harmonic expansion of the anisotropy in equatorial coordinates this corresponds to a filter of m = 0 multipole moments, which introduces a systematic underestimation of the observed anisotropy. In particular, the dipole anisotropy, which is a crucial observable for the study of cosmic ray diffusion in our Galaxy, can only be observed as a projection onto the celestial equator.
This has important consequences for the interpretation of experimental dipole data. If the dipole orientation changes with rigidity the projected dipole can exhibit rigidity modulations introduced by the projection on top of a true rigidity dependence of the dipole amplitude. If the dipole vector aligns with the celestial poles these projection effects can lead to a drastic reduction of the observed dipole accompanied by a phase-flip.
In addition, the limited integrated field of view of observatories affects the reconstruction of the multipole moments. As as consequence, even the m = 0 filtered anisotropy cannot be unambiguously reconstructed. This can introduce a large systematic uncertainty, in particular, for the low-multipole moments of the anisotropy. We have discussed strategies how to account for this effect in the multipole reconstruction.
We would like to conclude with a few remarks. In our anisotropy ansatz (2), we have not accounted for varia-tions in the primary cosmic-ray flux during the observation period. The solar potential, with a variation on time scales of 11 years, only affects low-energy data and can be neglected for cosmic rays in the TeV to PeV energy range (Munakata et al. 2010). On the other hand, the relative motion of the Earth in the solar system is expected to be visible as a solar dipole (Compton & Getting 1935). This effect would be visible if the data is binned in terms of solar time. Indeed, a solar dipole anisotropy has been observed at the level of 10 −4 in multi-TeV cosmic-ray data (Amenomori et al. 2004(Amenomori et al. , 2006Abdo et al. 2009;Abbasi et al. 2011Abbasi et al. , 2012Bartoli et al. 2015). However, the solar dipole is expected to average out for the data binned in sidereal time as long as the observation period covers an integer number of full years.
The method presented in this paper is designed for midlatitude observatories that are exposed to different parts of the celestial sky as the Earth rotates. The method is also applicable to the case of IceCube/IceTop located at the South Pole. However, for these observatories, the instantaneous field of view is identical to the time-integrated one. As a result, the iteration does not improve the estimates from direct integration or time scrambling. On the other hand, the special location in combination with a limited field of view makes these observatories particularly insensitive to low-multipoles due to projection effects of the anisotropy discussed in the text.
In this paper, we have only considered the case of the data analysis of individual observatories. A great advantage of the likelihood-based reconstruction method is the straightforward generalization to combined anisotropy studies of data sets from multiple observatories with overlapping integrated field of view, see e.g. (Díaz-Vélez, J. C. 2015). The expectation value (4) can then be simply generalized to a sum over data sets with individual detector exposures but same anisotropy, as long as the rigidity distributions of the data sets are very similar.
We would also like to emphasize that cosmic-ray observations via satellites like Fermi could have an advantage since the observatory can be tilted during the observation. In this case, it is possible to break the degeneracy between local acceptance effects and cosmicray anisotropy. In principle, it is then possible to have a full reconstruction of the anisotropy without projection effects, provided that additional systematic effects are under control.
The partition of the field of view into multiple sectors is at this stage arbitrary and should in general be guided by the property of the data. For instance, the Equal-Zenith-Angle Method divides the sky into ring-like sectors with