LARGE-SCALE DISTRIBUTION OF ARRIVAL DIRECTIONS OF COSMIC RAYS DETECTED ABOVE 10¹⁸ eV AT THE PIERRE AUGER OBSERVATORY

The energy estimator of the showers recorded by the SD array is provided by the signal at 1000 m from the shower core, S(1000). For any fixed energy, since the development of extensive air showers depends on the atmospheric pressure P and air density ρ, the corresponding S(1000) is sensitive to variations in pressure and air density. Systematic variations with time of S(1000) induce variations of the event rate that may distort the real dependence of the cosmic-ray intensity with right ascension. To cope with this experimental effect, the observed shower size S(1000), measured at the actual density ρ and pressure P, is related to the one S_atm(1000) that would have been measured at reference values ρ₀ and P₀ (The Pierre Auger Collaboration 2009):

$$\begin{eqnarray} S_{{\rm atm}}(1000)&=&[1-\alpha _P(\theta)(P-P_0)-\alpha _\rho (\theta)(\rho _d-\rho _0)\nonumber\\ && -\beta _\rho (\theta)(\rho -\rho _d)]S(1000). \end{eqnarray} \tag{ 1 }$$

The reference values are chosen as the average values at Malargüe (i.e., ρ₀ = 1.06 kg m⁻³ and P₀ = 862 hPa). ρ_d denotes here the average daily density at the time the event was recorded. The measured coefficients α_ρ, β_ρ and α_P—given in Table 1—give the influence on the shower sizes of the air density (and thus temperature) at long and short timescales on the Molière radius (and hence the lateral profiles of the showers) and of the pressure on the longitudinal development of air showers, respectively.

Table 1. Coefficients α_ρ, β_ρ, and α_P Used to Correct Shower Sizes for Atmospheric Effects on Shower Development, in Bins of $\sec {\theta }$

$\sec {\theta }$	αρ(kg−1 m3)	βρ(kg−1 m3)	αP(h Pa−1)
[1.0–1.2]	−9.7 10−1	−2.6 10−1	−4.4 10−4
[1.2–1.4]	−7.2 10−1	−2.2 10−1	−1.6 10−3
[1.4–1.6]	−5.4 10−1	−2.0 10−1	−2.3 10−3
[1.6–1.8]	−4.0 10−1	−4.3 10−2	−1.9 10−3
[1.8–2.0]	−1.5 10−1	−2.3 10−2	−2.8 10−3

Note. From The Pierre Auger Collaboration (2009).

Download table as:  ASCIITypeset image

Applying these corrections to the energy assignments of showers allows us to cancel spurious variations of the event rate in right ascension, where typical amplitudes amount to a few per thousand when considering data sets collected over full years.

The trajectories of charged particles in extensive air showers are curved in Earth's magnetic field, resulting in a broadening of the spatial distribution of particles in the direction of the Lorentz force. As the strength of the geomagnetic field component perpendicular to any arrival direction depends on both the zenith and azimuthal angles, the small changes of the density of particles at ground induced by the field break the circular symmetry of the lateral spread of the particles and thus induce a dependence of the shower size S(1000) at a fixed energy in terms of the azimuthal angle. Due to the steepness of the energy spectrum, such an azimuthal dependence translates into azimuthal modulations of the estimated cosmic-ray event rate at a given S(1000). To eliminate these effects, the observed shower size S(1000) is related to the one that would have been observed in the absence of geomagnetic field S_geom(1000) (The Pierre Auger Collaboration 2011b):

$$\begin{eqnarray} &&S_{{\rm geom}}(1000)=\left[1-g_1\cos ^{-g_2}{(\theta)}\sin ^2{(\widehat{\textbf {u},\textbf {b}})}\right]S(1000), \end{eqnarray} \tag{ 2 }$$

where g₁ = (4.2 ± 1)10⁻³, g₂ = 2.8 ± 0.3, and $\textbf {u}$ and $\textbf {b}=\textbf {B}/\Vert \textbf {B} \Vert$ denote the unit vectors in the shower direction and the geomagnetic field direction, respectively. At a zenith angle θ = 55°, the amplitude of the asymmetry in azimuth already amounts to ≃ 2%, which is why we restrict the present analysis to zenith angles smaller than this value. Carrying out these corrections is thus critical for performing large-scale anisotropy measurements in declination.

Once the influence on S(1000) of weather and geomagnetic effects are accounted for, the dependence of S(1000) on zenith angle, due to the attenuation of the shower and geometrical effects, is extracted from the data using the constant intensity cut method (The Pierre Auger Collaboration 2008). The attenuation curve CIC(θ) is fitted with a second-order polynomial in x = cos ²(θ) − cos ²(38°):CIC(θ) = 1 + ax + bx². The angle 38° is chosen as a reference to convert S(1000) to $S_{38^\circ }=S(1000)/{\it CIC}(\theta)$. $S_{38^\circ }$ may be regarded as the signal that would have been expected had the shower arrived at 38°. The values of the parameters a = 0.94 ± 0.03 and b = −0.95 ± 0.05 are deduced for $S_{38^\circ }=22$VEM,¹⁰² which corresponds to an energy of about 4 EeV—just above the threshold energy for full efficiency. The differences between these parameters and previous reports will be discussed in Section 6.

Finally, the sub-sample of events recorded by both the fluorescence telescopes and the SD array is used to establish the relationship between the energy reconstructed with the fluorescence telescopes E_FD and $S_{38^\circ }$:E_FD = AS^B_38°. The resulting parameters from the data fit are A = (1.68 ± 0.05) × 10⁻¹ EeV and B = 1.030 ± 0.009, in good agreement with the recent report given in Pesce et al. (2011). The energy scale inferred from this data sample is applied to all showers detected by the SD array.

The choice of the fiducial cut to select high-quality events allows the precise determination of the geometric directional aperture per cell as a_cell(θ) = 1.95cos θ km² (The Pierre Auger Collaboration 2010b). It also allows us to exploit the regularity of the array for obtaining its geometric directional aperture as a simple multiple of a_cell(θ) (The Pierre Auger Collaboration 2010b). The number of elemental cells n_cell(t) is accurately monitored every second at the Observatory. To search for celestial large-scale anisotropies, it is mandatory to account for the modulation imprinted by the variations of n_cell(t) in the expected number of events at the sidereal periodicity T_sid. Within each sidereal day, and in the same way as in The Pierre Auger Collaboration (2011a), we denote by α⁰ the local sidereal time and express it in hours or in radians, as appropriate. For practical reasons, α⁰ is chosen so that it is always equal to the right ascension of the zenith at the center of the array. As a function of α⁰, the total number of elemental cells N_cell(α⁰) and its associated relative variations ΔN_cell(α⁰) are then obtained from

$$\begin{equation} N_{\mathrm{cell}}(\alpha ^0)=\sum _{j}n_{\mathrm{cell}}(\alpha ^0+jT_{{\rm sid}}), \quad \Delta N_{\mathrm{cell}}(\alpha ^0)=\frac{N_{\mathrm{cell}}(\alpha ^0)}{\left<N_{\mathrm{cell}}\right>_{\alpha ^0}}, \end{equation} \tag{ 3 }$$

with $\left<N_{\mathrm{cell}}\right>_{\alpha ^0}=1/T_{{\rm sid}}\int _0^{T_{{\rm sid}}}{d}\alpha ^0N_{\mathrm{cell}}(\alpha ^0)$. In the same way as in The Pierre Auger Collaboration (2011a), the small modulation of the expected number of events in right ascension induced by those variations will be accounted for by weighting each event k with a factor inversely proportional to ΔN_cell(α⁰_k) when estimating the anisotropy parameters in Section 5. Placing such time dependences in the event weights allows us to remove the modulations in time imprinted by the growth of the array and the dead times for each detector.

At any time, the effective directional aperture of the SD array is controlled by the geometric one and by the detection efficiency function (θ, φ, E). For each elemental cell, the directional exposure in celestial coordinates is then simply obtained through the integration over local sidereal time of x⁽ⁱ⁾(α⁰) × a_cell(θ) × (θ, φ, E), where x⁽ⁱ⁾(α⁰) is the operational time of the cell (i). Actually, since the small modulations in time imprinted in the event counting rate by experimental effects will be accounted for by means of the weighting procedure just described when searching for anisotropies, the small variations in local sidereal time for each x⁽ⁱ⁾(α⁰) can be neglected in calculating ω. The zenith and azimuth angles are related to the declination and the right ascension through

$$\begin{eqnarray} \cos {\theta }&=&\sin {\delta }\sin {\ell _\mathrm{site}}+\cos {\delta }\cos {\ell _\mathrm{site}}\cos {(\alpha -\alpha ^0)},\nonumber \\ \tan {\varphi }&=&\frac{\cos {\delta }\sin {\ell _\mathrm{site}}\cos {(\alpha -\alpha ^0)}-\sin {\delta }\cos {\ell _\mathrm{site}}}{\cos {\delta }\sin {(\alpha -\alpha ^0)}}, \end{eqnarray} \tag{ 4 }$$

with ℓ_site the mean latitude of the Observatory. Since both θ and φ depend only on the difference α − α⁰, the integration over α⁰ can then be substituted for an integration over the hour angle α' = α − α⁰ so that the directional exposure actually does not depend on right ascension when the x⁽ⁱ⁾ are assumed local sidereal time independent:

$$\begin{eqnarray} \omega (\delta,E) &=& \sum _{i=1}^{n_{\mathrm{cell}}}x^{(i)}\int _0^{24h} {d}\alpha ^\prime \,a_{\mathrm{cell}}{(\theta (\alpha ^\prime,\delta))} \nonumber\\ && \times \epsilon (\theta (\alpha ^\prime,\delta),\varphi (\alpha ^\prime,\delta),E). \end{eqnarray} \tag{ 5 }$$

Above 3 EeV, this integration can be performed analytically (Sommers 2001). Below 3 EeV, the non-saturation of the detection efficiency makes the directional exposure lower. The next sections are dedicated to the determination of (θ, φ, E).

To determine the detection efficiency function, a natural method would be to generate showers by means of Monte Carlo simulations and to calculate the ratio of the number of triggered events to the total simulated. However, there are discrepancies in the predictions of the hadronic interaction model regarding the number of muons in shower simulations and what is found in our data (Engel et al. 2007). This prevents us from relying on this method for obtaining the detection efficiency to the required accuracy.

We adopt here instead an empirical approach, based on the quasi-invariance of the zenithal distribution to large-scale anisotropies for zenith angles less than ≃ 60° and for any Observatory whose latitude is far from the poles of the Earth. For full efficiency, the distribution in zenith angles dN/dθ is proportional to sin θcos θ for solid angle and geometry reasons, so that the distribution in dN/dsin ²θ is uniform. Consequently, below full efficiency, any significant deviation from a uniform behavior in the dN/dsin ²θ distribution provides an empirical measurement of the zenithal dependence of the detection efficiency. The quasi-invariance of dN/dsin ²θ to large-scale anisotropies is demonstrated in Appendix A.

Based on this quasi-invariance, the detection efficiency averaged over the azimuth can be estimated from

$$\begin{equation} \left<\epsilon (\theta,\varphi,E)\right>_{\varphi } = \frac{1}{\mathcal {N}}\frac{{d}N(\sin ^2{\theta },E)}{{d}\sin ^2{\theta }}, \end{equation} \tag{ 6 }$$

where the notation <· > _φ stands for the average over φ and the constant $\mathcal {N}$ is the number of events that would have been observed at energy E and for any sin ²θ value in case of full efficiency for an energy spectrum dN/dE = 40(E/EeV)^−3.27 km⁻² yr⁻¹ sr⁻¹ EeV⁻¹—as measured between 1 and 4 EeV (The Pierre Auger Collaboration 2010a). Consequently, for each zenith angle, this empirical measurement of the efficiency provides an estimate relative to the overall spectrum of cosmic rays. In particular, since it is applied to all events detected at energy E without distinction based on the primary mass of cosmic rays, this technique does not provide the mass dependence of the detection efficiency. For that reason, the anisotropy searches reported in Section 5 pertain to the whole population of cosmic rays, whether this population consists of a single primary mass or a mixture of several elements.

Results are shown in Figure 1 for four different energies.¹⁰⁴ At 4 EeV, a uniform behavior around 1 is observed, though it is quite noisy due to the reduced statistics. This uniform behavior is consistent with full efficiency at this energy, as expected. Note that some values are greater than 1 for energies close to or higher than 3 EeV because of the empirical method of measuring the efficiency relative to the overall spectrum of cosmic rays. At 2 EeV, a loss of efficiency is observed for vertical showers due to the attenuation of the electromagnetic component of the showers. Up to ≃ 40°, the detection efficiency steadily increases because the projected area of showers at ground gets larger with zenith angle. Above ≃ 40°, the rapid increase of the slant depth then makes the attenuation of the electromagnetic component stronger, but the muonic component of showers becomes dominant and ensures a high detection efficiency. At lower energies, the number of muons is, in contrast, too low to significantly impact the detection efficiency above ≃ 40°–45°, so that a clear decrease is observed at high zenith angles. In the following, we use parameterizations obtained by fitting each distribution with a fourth-order polynomial function in sin ²θ, which is sufficient to reproduce the main details as illustrated in Figure 1.

Zoom In Zoom Out Reset image size

In addition to the effects on the energy determination presented in Section 3.2, geomagnetic effects also affect the detection efficiency for showers with energies below 3 EeV. This is because under any incident angles (θ, φ), a shower with an energy E triggers the SD array with a probability associated with its size which is a function of azimuth because of the geomagnetic effects¹⁰⁵: E × (1 + Δ(θ, φ))^B. Above 1 EeV, this effect is in fact the main source of azimuthal dependence of the detection efficiency, so that to first order in Δ(θ, φ), (θ, φ, E) can be estimated as

$$\begin{eqnarray} \epsilon (\theta,\varphi,E) &=& \frac{1}{\mathcal {N}}\frac{{d}N(\sin ^2{\theta },E(1+\Delta (\theta,\varphi))^B)}{{d}\sin ^2{\theta }} \nonumber \\ &\simeq &\left<\epsilon (\theta,\varphi,E)\right>_{\varphi }+\frac{BE\Delta (\theta,\varphi)}{\mathcal {N}}\frac{\partial \left<\epsilon (\theta,\varphi,E)\right>_{\varphi }}{\partial E}.\quad\quad \end{eqnarray} \tag{ 7 }$$

The correction to the detection efficiency induced by geomagnetic effects, and, in particular, the azimuthal dependence, is thus straightforward to implement from the knowledge of <(θ, φ, E) > _φ. An example of such an azimuthal dependence is shown in the left panel of Figure 2, for E = 1 EeV and θ = 55°. The modulation reflects the one due to the energy determination: the detection efficiency is lowered in the directions where the uncorrected energies are underestimated due to geomagnetic effects, and the efficiency is higher where energies are overestimated. The maximal contrast of such azimuthal modulations is displayed in the right panel as a function of the zenith angle, for three different energies. At 2 EeV, the amplitude slightly increases up to ≃ 35°, staying below ≃ 0.1%, and then decreases and even cancels due to the saturation of the detection efficiency. In contrast, when decreasing in energy, the relative amplitude largely increases with the zenith angle due to the increase of the derivative term, reaching ≃ 1.7% for θ = 55° and E = 1 EeV.

Zoom In Zoom Out Reset image size

The altitudes above sea level of the water-Cherenkov detectors are displayed in Figure 3 with color coding. The coordinates are in a Cartesian system whose origin is defined at the "center" of the Observatory site. The Andes ridge building up in the western and northwestern direction can be seen. A slightly tilted SD array gives rise to a small azimuthal asymmetry, and consequently slightly modifies the directional exposure with respect to Equation (5) through small changes of the geometric directional aperture. This modification is twofold: the tilt changes the geometric factor (cos θ) of the projected surface under incidence angles (θ, φ) and also induces a compensating effect below full efficiency by slightly varying the detection efficiency with the azimuth angle φ.

Zoom In Zoom Out Reset image size

Denoting $\mathbf {n_\perp ^{(i)}}$ the normal vector to each elemental cell, the geometric directional aperture per cell is no longer simply given by cos θ but now depends on both θ and φ:

$$\begin{eqnarray} a_{\mathrm{cell}}^{(i)}(\theta,\varphi)&=&1.95\, \mathbf {n}\cdot \mathbf {n_\perp ^{(i)}} \simeq 1.95 \nonumber\\ && \times \big[1+\zeta ^{(i)}\tan {\theta }\cos {\big(\varphi -\varphi _0^{(i)} \big)} \big] \cos {\theta }, \end{eqnarray} \tag{ 8 }$$

where ζ⁽ⁱ⁾ and φ⁽ⁱ⁾₀ are the zenith and azimuth angles of $\mathbf {n_\perp ^{(i)}}$. It is actually this latter expression a_cell which has to be inserted into Equation (5) to calculate the directional exposure. Overall, the average tilt of the SD array is ζ^eff ≃ 02, and induces a dipolar asymmetry in azimuth with a maximum in the downhill direction φ^eff₀ ≃ 0° and with an amplitude increasing with the zenith angle as ≃ 0.3%tan θ.

Below 3 EeV, the tilt of the array induces an additional variation of the detection efficiency with azimuth. This is because the effective separation between detectors for a given zenith angle now depends on the azimuth. Since, for a given zenith angle, the SD array seen by showers coming from the uphill direction is denser than that for those coming from the downhill direction, the detection efficiency is higher in the uphill direction. Parameterizing the energy dependence of as E³/(E³ + E³_0.5), we show in Appendix B that the change in the detection efficiency can be estimated as

$$\begin{equation} \Delta \epsilon _{\mathrm{tilt}}(\theta,\varphi,E) = \frac{E^3\big(E_{0.5}^3-{E_{0.5}^{\mathrm{tilt}}}^3(\theta,\varphi)\big)}{\big(E^3+E_{0.5}^3\big)\big(E^3+{E_{0.5}^{\mathrm{tilt}}}^3(\theta,\varphi)\big)}, \end{equation} \tag{ 9 }$$

where E^tilt_0.5(θ, φ) is related to E_0.5 through

$$\begin{equation} E_{0.5}^{\mathrm{tilt}}(\theta,\varphi)\simeq E_{0.5}\times \big[1+\zeta ^{\mathrm{eff}}\tan {\theta }\cos {\big(\varphi -\varphi _0^{\mathrm{eff}} \big)}\big]^{3/2}. \end{equation} \tag{ 10 }$$

Around 1 EeV, this correction tends to compensate the pure geometrical effect described above, and even overcompensates it at lower energies.

This spatial extension of the SD array is such that the range of latitudes covered by all cells reaches ≃ 05. This induces a slightly different directional exposure between the cells located at the northern part of the array and the ones located at the southern part. This spatial extension can be accounted for to calculate the overall directional exposure using the cell latitudes ℓ⁽ⁱ⁾_cell instead of the mean site one in the transformations from local to celestial angles in Equation (4).

In the same way as geomagnetic effects, weather effects can also affect the detection efficiency for showers with energies below 3 EeV. However, above 1 EeV, we have shown in The Pierre Auger Collaboration (2011a) that as long as the analysis covers an integer number of years with almost equal exposure in every season, the amplitude of the spurious modulation in right ascension induced by this effect is small enough to be neglected when performing anisotropy analyses at the present level of sensitivity.

Accounting for all effects, the final expression to calculate the directional exposure is slightly modified with respect to Equation (5):

$$\begin{eqnarray} \omega (\delta,E) &=& \sum _{i=1}^{n_{\mathrm{cell}}}x^{(i)}\int _0^{24h} {d}\alpha ^\prime \,a_{\mathrm{cell}}^{(i)}{(\theta,\varphi)}\nonumber\\ && \times [\epsilon (\theta,\varphi,E)+\Delta \epsilon _{\mathrm{tilt}}(\theta,\varphi,E)], \end{eqnarray} \tag{ 11 }$$

where both θ and φ depend on α', δ, and ℓ⁽ⁱ⁾_cell. The resulting dependence on declination is displayed in Figure 4 for three different energies. Down to 1 EeV, the detection efficiency at high zenith angles is high enough that the equatorial south pole is visible at any time and hence constitutes the direction of maximum of exposure. For a wide range of declinations between ≃ − 89° and ≃ − 20°, the directional exposure is ≃ 2500 km² yr at 1 EeV, and ≃ 3500 km² yr for any energy above full efficiency. Then, at higher declinations, it smoothly falls to zero, with no exposure above ≃ 20° declination.

Zoom In Zoom Out Reset image size

The average expected number of events within any solid angle and any energy range can be recovered by integrating the directional exposure over the solid angle considered and the cosmic-ray energy spectrum in the corresponding energy range. Note that the rapid variation of the exposure close to the South Pole on an angular scale of the order of the angular resolution has no influence on the event counting rate, due to the quasi-zero solid angle in that particular direction. Consequently, though the exposure around the South Pole could be affected by small changes of the detection efficiency around θ = 55°, the results presented in next sections are on the other hand not affected by the exact value of the exposure for declinations a few degrees away from the South Pole.

Any angular distribution over the sphere $\Phi (\mathbf {n})$ can be decomposed in terms of a multipolar expansion:

$$\begin{equation} \Phi (\mathbf {n})=\sum _{\ell \ge 0}\sum _{m=-\ell }^{\ell } a_{\ell m}Y_{\ell m}(\mathbf {n}), \end{equation} \tag{ 12 }$$

where $\mathbf {n}$ denotes a unit vector taken in equatorial coordinates. The customary recipe to extract each multipolar coefficient makes use of the completeness relation of spherical harmonics:

$$\begin{equation} a_{\ell m}=\int _{4\pi } {d}\Omega \Phi (\mathbf {n})Y_{\ell m}(\mathbf {n}), \end{equation} \tag{ 13 }$$

where the integration is over the entire sphere of directions $\mathbf {n}$. Any anisotropy fingerprint is encoded in the a_ℓm spherical harmonic coefficients. Variations on an angular scale of Θ radians contribute amplitude in the ℓ ≃ 1/Θ modes.

However, in the case of partial sky coverage, the solid angle in the sky where the exposure is zero makes it impossible to estimate the multipolar coefficients a_ℓm in this way. This is because the unseen solid angle prevents one from making use of the completeness relation of the spherical harmonics (Sommers 2001). Since the observed arrival direction distribution is in this case the combination of the angular distribution $\Phi (\mathbf {n})$ and of the directional exposure function $\omega (\mathbf {n})$, the integration performed in Equation (13) does not allow any longer the extraction of the multipolar coefficients of $\Phi (\mathbf {n})$, but only the ones of $\omega (\mathbf {n}) \Phi (\mathbf {n})$ (Billoir & Deligny 2008)¹⁰⁶:

$$\begin{eqnarray} b_{\ell m}&=&\int _{\Delta \Omega } {d}\Omega \omega (\mathbf {n})\Phi (\mathbf {n})Y_{\ell m}(\mathbf {n}) \nonumber \\ &=&\sum _{\ell ^\prime \ge 0}\sum _{m^\prime =-\ell ^\prime }^{\ell ^\prime } a_{\ell ^\prime m^\prime }\int _{\Delta \Omega } {d}\Omega \omega (\mathbf {n})Y_{\ell ^\prime m^\prime }(\mathbf {n})Y_{\ell m}(\mathbf {n}). \end{eqnarray} \tag{ 14 }$$

Formally, the a_ℓm coefficients appear related to the b_ℓm ones through a convolution such that $b_{\ell m}=\sum _{\ell ^\prime \ge 0}\sum _{m^\prime =-\ell ^\prime }^{\ell ^\prime }[K]_{\ell m}^{\ell ^\prime m^\prime } a_{\ell ^\prime m^\prime }$. The matrix K, which imprints the interferences between modes induced by the non-uniform and partial coverage of the sky, is entirely determined by the directional exposure. The relationship established in Equation (14) is valid for any exposure function $\omega (\mathbf {n})$.

Meanwhile, the observed arrival direction distribution, $\overline{{d}N}(\mathbf {n})/{d}\Omega$, provides a direct estimation of the b_ℓm coefficients through (hereafter, we use an overline to indicate the estimator of any quantity)

$$\begin{equation} \overline{b}_{\ell m}=\int _{\Delta \Omega } {d}\Omega \frac{\overline{{d}N}(\mathbf {n})}{{d}\Omega } Y_{\ell m}(\mathbf {n}), \end{equation} \tag{ 15 }$$

where the distribution $\overline{{d}N}(\mathbf {n})/{d}\Omega$ of any set of N arrival directions $\lbrace \mathbf {n}_1,..., \mathbf {n}_N\rbrace$ can be modeled as a sum of Dirac functions on the sphere. Then, if the multipolar expansion of the angular distribution $\Phi (\mathbf {n})$ is bounded to ℓ_max, that is, if the $\Phi (\mathbf {n})$ has no higher moments than ℓ_max, the first b_ℓm coefficients with ℓ ⩽ ℓ_max are related to the non-vanishing a_ℓm by the square matrix $K_{\ell _{\mathrm{max}}}$ truncated to ℓ_max. Inverting this truncated matrix allows us to recover the underlying a_ℓm from the measured b_ℓm (with ℓ ⩽ ℓ_max):

$$\begin{equation} \overline{a}_{\ell m}=\sum _{\ell ^\prime =0}^{\ell _{\mathrm{max}}}\sum _{m^\prime =-\ell ^\prime }^{\ell ^\prime } \big[K^{-1}_{\ell _{\mathrm{max}}}\big]_{\ell m}^{\ell ^\prime m^\prime } \overline{b}_{\ell ^\prime m^\prime }. \end{equation} \tag{ 16 }$$

In the case of small anisotropies (|a_ℓm|/a₀₀ ≪ 1), the resolution on each recovered $\overline{a}_{\ell m}$ coefficient is proportional to $([K^{-1}_{\ell _{\mathrm{max}}}]_{\ell m}^{\ell m})^ {0.5}$ (Billoir & Deligny 2008):

$$\begin{equation} \sigma _{\ell m}=\Big (\big[K^{-1}_{\ell _{\mathrm{max}}}\big]_{\ell m}^{\ell m} \overline{a}_{00}\Big)^{0.5}. \end{equation} \tag{ 17 }$$

The dependence on ℓ_max of the coefficients of $K^{-1}_{\ell _{\mathrm{max}}}$ induces an intrinsic indeterminacy of each recovered coefficient $\overline{a}_{\ell m}$ as ℓ_max is increasing. This is nothing else but the mathematical translation of it being impossible to know the angular distribution of cosmic rays in the uncovered region of the sky.

Henceforth, we adapt this general formalism to the search for anisotropies in Auger data in different energy intervals. We assume that the energy dependence of the angular distribution of cosmic rays is smooth enough that the multipolar coefficients can be considered constant for any energy E within a narrow interval ΔE. The directional exposure is hereafter considered as independent of the right-ascension, as defined in Section 4. Within an energy interval ΔE, the expected arrival direction distribution thus reads

$$\begin{equation} \frac{{d}N(\mathbf {n})}{{d}\Omega }\propto \tilde{\omega }(\delta) \sum _{\ell \ge 0}\sum _{m=-\ell }^{\ell } a_{\ell m}Y_{\ell m}(\mathbf {n}), \end{equation} \tag{ 18 }$$

where $\tilde{\omega }(\delta)$ is the effective directional exposure for the energy interval ΔE. For convenience, this latter function is normalized such that

$$\begin{equation} \tilde{\omega }(\delta)= \frac{\displaystyle \int _{\Delta E} {d}E E^{-\gamma }\omega (\delta,E)}{\displaystyle \max _\delta \bigg [\int _{\Delta E} {d}E E^{-\gamma }\omega (\delta,E)\bigg ]}, \end{equation} \tag{ 19 }$$

with γ the spectral index in the considered energy range. This dimensionless function provides, for any direction on the sky, the effective directional exposure in the energy range ΔE at that direction, relative to the largest directional exposure on the sky. This is actually the relevant quantity which enters into Equation (14) for the analyses presented below. Note that for a directional exposure independent of the right ascension, the coefficients [K]^ℓ'm'_ℓm are proportional to δ^m'_m, i.e., different values of m are not mixed in the matrix. The observed arrival direction distribution, $\overline{{d}N}(\mathbf {n})/{d}\Omega$, is here modeled as a sum of Dirac functions on the sphere weighted by the factor ΔN⁻¹_cell(α⁰_k) for each event recorded at local sidereal time α⁰_k, as described in Section 4.1 to correct for the slightly non-uniform directional exposure in right ascension. In this way, the integration in Equation (14) yields to

$$\begin{equation} \overline{b}_{\ell m}=\sum _{k=1}^{N} \frac{Y_{\ell m}(\mathbf {n}_k)}{\Delta N_{\mathrm{cell}}\big(\alpha ^0_{k}\big)}. \end{equation} \tag{ 20 }$$

The multipolar coefficients $\overline{a}_{\ell m}$ are then recovered by means of Equation (16). Given the exposure functions described in Section 4, the resolution on each recovered coefficient, encoded in Equation (17), is degraded by a factor larger than 2 each time ℓ_max is increased by 1. This prevents the recovery of each coefficient with good accuracy as soon as ℓ_max ⩾ 3, since, for ℓ_max = 3 for instance, our current statistics would only allow us to probe dipole amplitudes at the 10% level. Consequently, in the following, we restrict ourselves to reporting results on individual coefficients obtained when assuming a dipolar distribution (ℓ_max = 1) and a quadrupolar distribution (ℓ_max = 2). Meanwhile, due to the interference between modes induced by the non-uniform and partial sky coverage, it is important to stress again that each multipolar coefficient recovered under the assumption of a particular bound ℓ_max might be biased if the underlying angular distribution of cosmic rays is not bounded to ℓ_max. Given the directional exposure functions considered in this study, this effect can be important only if the angular distribution has in fact significant moments of order ℓ_max + 1.

As outlined in the Introduction, a measurable dipole is regarded as a likely possibility in many scenarios for the origin of cosmic rays at EeV energies. Assuming that the angular distribution of cosmic rays is modulated by a pure dipole, the intensity $\Phi (\mathbf {n})$ can be parameterized in any direction $\mathbf {n}$ as

$$\begin{equation} \Phi (\mathbf {n})=\frac{\Phi _0}{4\pi } (1+r \mathbf {d}\cdot \mathbf {n}), \end{equation} \tag{ 21 }$$

where $\mathbf {d}$ denotes the dipole unit vector. The dipole pattern is here fully characterized by a declination δ_d, a right ascension α_d, and an amplitude r corresponding to the maximal anisotropy contrast:

$$\begin{equation} r=\frac{\Phi _{\mathrm{max}}-\Phi _{\mathrm{min}}}{\Phi _{\mathrm{max}}+\Phi _{\mathrm{min}}}. \end{equation} \tag{ 22 }$$

The estimation of these three coefficients is straightforward from the estimated spherical harmonic coefficients $\overline{a}_{1m}$: $\overline{r}=[3(\overline{a}_{10}^2+\overline{a}_{11}^2+\overline{a}_{1-1}^2)]^{0.5}/\overline{a}_{00}$, $\overline{\delta }=\arcsin {(\sqrt{3}\overline{a}_{10}/\overline{a}_{00}\overline{r})}$, and $\overline{\alpha }=\arctan {(\overline{a}_{1-1}/\overline{a}_{11})}$. Uncertainties on $\overline{r}$, $\overline{\delta }$, and $\overline{\alpha }$ are obtained from the propagation of uncertainties on each recovered $\overline{a}_{1m}$ coefficient (cf. Equation (17)). Under an underlying isotropic distribution, and for an axisymmetric directional exposure around the axis defined by the North and South equatorial poles, the probability density function of $\overline{r}$ is given by (The Pierre Auger Collaboration 2011b):

$$\begin{equation} p_R(\overline{r})=\frac{\overline{r}}{\sigma \sqrt{\sigma _z^2-\sigma ^2}} \mathrm{erfi}\bigg (\frac{\sqrt{\sigma _z^2-\sigma ^2}}{\sigma \sigma _z}\frac{\overline{r}}{\sqrt{2}}\bigg) \exp {\bigg (-\frac{\overline{r}^2}{2\sigma ^2}\bigg)}, \vspace*{4pt} \end{equation} \tag{ 23 }$$

where erfi(z) = erf(iz)/i, $\sigma =\sqrt{3}\sigma _{11}/\overline{a}_{00}$, and $\sigma _z=\sqrt{3}\sigma _{10}/\overline{a}_{00}$. The probability $P_R(>\overline{r})$ that an amplitude equal to or larger than $\overline{r}$ arises from a statistical fluctuation of an isotropic distribution is then obtained by integrating p_R above $\overline{r}$:

$$\begin{equation} P_R(>\overline{r})=\mathrm{erfc}\bigg (\!\frac{\overline{r}}{\sqrt{2}\sigma _z} \!\bigg)+\mathrm{erfi}\bigg (\!\frac{\sqrt{\sigma _z^2-\sigma ^2}}{\sigma \sigma _z}\frac{\overline{r}}{\sqrt{2}} \!\bigg) \exp {\bigg (\!-\frac{\overline{r}^2}{2\sigma ^2} \!\bigg)}. \vspace*{4pt} \end{equation} \tag{ 24 }$$

The reconstructed amplitudes $\overline{r}(E)$ and corresponding directions are shown in Figure 5 with the associated uncertainties, as a function of the energy. The directions are drawn in azimuthal projection, with the equatorial South Pole located at the center and the right ascension going from 0 to 360° clockwise. In the left panel, the 99% CL upper bounds on the amplitudes that would result from fluctuations of an isotropic distribution are indicated by the dotted line (i.e., the amplitudes $\overline{r}_{99}(E)$ such that $P_R(>\overline{r}_{99}(E))=0.01$). One can see that within the statistical uncertainties, there is no strong evidence of any significant signal.

Zoom In Zoom Out Reset image size

The reconstructed declinations $\overline{\delta }$ and right ascensions $\overline{\alpha }$ are shown separately in Figure 6. Both quantities are expected to be randomly distributed in case of independent samples whose parent distribution is isotropic. In our previous report on the first harmonic analysis in right ascension (The Pierre Auger Collaboration 2011a), we pointed out the intriguing smooth alignment of the phases in right ascension as a function of the energy, and noted that such a consistency of phases in adjacent energy intervals is expected to manifest with a smaller number of events than those required for the detection of amplitudes standing out significantly above the background noise in the case of a real underlying anisotropy. This motivated us to design a prescription aimed at establishing at 99% CL whether this consistency in phases is real, using the exact same analysis as the one reported in The Pierre Auger Collaboration (2011a). The prescribed test will end once the total exposure since 2011 June 25 reaches 21,000 km² yr sr. The smooth fit to the data of The Pierre Auger Collaboration (2011a) is shown as a dashed line in the right panel of Figure 6, restricted to the energy range considered here. Though the phase between 4 and 8 EeV is poorly determined due to the corresponding direction in declination pointing close to the equatorial South Pole, it is noteworthy that a consistent smooth behavior is observed using the analysis presented here and applied to a data set containing two additional years of data. It is also interesting to see in the left panel that all reconstructed declinations are in the equatorial Southern Hemisphere.

Zoom In Zoom Out Reset image size

For completeness, significance sky maps are displayed in Figure 7 in equatorial coordinates and using a Mollweide projection, for the four energy ranges. The Galactic plane and Galactic center are also depicted as the dotted line and the star. Significances are calculated using the Li and Ma estimator (Li & Ma 1983). This widely used estimator of significance, S, properly accounts for the fluctuations of the background and of an eventual signal in any angular region searched.¹⁰⁷ If no signal is present, the variable S is nearly normally distributed even for small count numbers, so that positive values of S can be interpreted as the number of standard deviations of any excess in the sky. Also, for negative values of S, −S can be interpreted as the number of standard deviations of any deficit in the sky. The maps show the overdensities obtained in circular windows of radius Θ = 1 radian, to better exhibit possible dipolar-like structures. The directions of the reconstructed dipoles are also shown, with their associated uncertainties (thick circles).

Zoom In Zoom Out Reset image size

Finally, since some consistency is observed both in declination and right ascension as a function of energy, the use of larger energy intervals and/or energy thresholds may help to pick up a significant signal above the background level. The amplitudes of the dipole are shown in Figure 8 for two energy intervals (1 < E/[EeV] < 4 and E > 4 EeV) and as a function of energy thresholds. This does not provide any further evidence for significant anisotropies.

Zoom In Zoom Out Reset image size

Any excesses along a plane would show up as a prominent quadrupole moment. Such excesses are plausible for instance at EeV energies in the case of an emission of light EeV cosmic rays from sources preferentially located in the Galactic disk, or at higher energies from sources preferentially located in the super-Galactic plane. Consequently, a measurable quadrupole may be regarded as an interesting outcome of an anisotropy search at ultrahigh energies.

Assuming now that the angular distribution of cosmic rays is modulated by a dipole and a quadrupole, the intensity $\Phi (\mathbf {n})$ can be parameterized in any direction $\mathbf {n}$ as

$$\begin{equation} \Phi (\mathbf {n})=\frac{\Phi _0}{4\pi } \bigg (1+r \mathbf {d}\cdot \mathbf {n} +\frac{1}{2}\sum _{i,j}Q_{ij}n_in_j \bigg), \end{equation} \tag{ 25 }$$

where $\mathbf {Q}$ is a traceless and symmetric second order tensor. Its five independent components are determined in a straightforward way from the ℓ = 2 spherical harmonic coefficients a_2m. Denoting by λ₊, λ₀, λ₋ the three eigenvalues of $\mathbf {Q}/2$ (λ₊ being the highest one and λ₋ the lowest one) and $\mathbf {q_+},\mathbf {q_0},\mathbf {q_-}$ the three corresponding unit eigenvectors, the intensity can be parameterized in a more intuitive way as

$$\begin{equation} \Phi (\mathbf {n})=\frac{\Phi _0}{4\pi } (1\,{+}\,r \mathbf {d}\cdot \mathbf {n} \,{+}\,\lambda _+(\mathbf {q_+}\cdot \mathbf {n})^2 \,{+}\,\lambda _0(\mathbf {q_0}\cdot \mathbf {n})^2 \,{+}\,\lambda _-(\mathbf {q_-}\cdot \mathbf {n})^2). \end{equation} \tag{ 26 }$$

It is then convenient to define the quadrupole amplitude β as

$$\begin{equation} \beta \equiv \frac{\lambda _+-\lambda _-}{2+\lambda _++\lambda _-}. \end{equation} \tag{ 27 }$$

In case of a pure quadrupolar distribution (i.e., in the absence of dipole), β is nothing else but the customary measure of maximal anisotropy contrast:

$$\begin{equation} r=0\Rightarrow \beta =\frac{\lambda _+-\lambda _-}{2+\lambda _++\lambda _-}=\frac{\Phi _{\mathrm{max}}-\Phi _{\mathrm{min}}}{\Phi _{\mathrm{max}}+\Phi _{\mathrm{min}}}. \end{equation} \tag{ 28 }$$

Hence, any quadrupolar pattern can be fully described by two amplitudes (β, λ₊) and three angles: (δ₊, α₊) which define the orientation of $\mathbf {q_+}$ and (α₋) which defines the direction of $\mathbf {q_-}$ in the orthogonal plane to $\mathbf {q_+}$. The third eigenvector $\mathbf {q_0}$ is orthogonal to $\mathbf {q_+}$ and $\mathbf {q_-}$, and its corresponding eigenvalue λ₀ is such that the traceless condition is satisfied: λ₊ + λ₋ + λ₀ = 0. Though the probability density functions of the estimated quadrupole amplitudes $(\overline{\beta },\overline{\lambda }_+)$ can be in principle calculated in the same way as in the case of the estimated dipole amplitude $(\overline{r})$, expressions are much more complicated to obtain even semi-analytically and we defer hereafter to Monte Carlo simulations to tabulate the distributions.

The amplitudes $\overline{r}(E)$, $\overline{\lambda }_+(E),$ and $\overline{\beta }(E)$ are shown in Figure 9 as functions of energy. Dipole amplitudes are compatible with expectations from isotropy. Compared to the results on the dipole obtained in previous section for ℓ_max = 1, the sensitivity is now degraded by a factor larger than 2 as expected from the dependence of the resolution σ_ℓm on ℓ_max (cf. Equation (17)). In the same way as for dipole amplitudes, the 99% CL upper bounds on the quadrupole amplitudes that could result from fluctuations of an isotropic distribution are indicated by the dashed lines. They correspond to the amplitudes $\overline{\lambda }_{+,99}(E)$ and $\overline{\beta }_{99}(E)$ such that the probabilities $P_{\Lambda _+}(>\overline{\lambda }_{+,99}(E))$ and $P_{B}(>\overline{\beta }_{99}(E))$ arising from statistical fluctuations of isotropy are equal to 0.01. Here, both distributions $P_{\Lambda _+}$ and P_B are sampled from Monte Carlo simulations. Throughout the energy scan, there is no evidence for anisotropy. The largest deviation from isotropic expectations occurs between 2 and 4 EeV, where the amplitude $\overline{\lambda }_+$ lies just above $\overline{\lambda }_{+99}$.

Zoom In Zoom Out Reset image size

Understanding the influence of the shower size corrections for geomagnetic effects is critical to get unbiased estimates of anisotropy parameters. Without accounting for these effects, an increase of the event rate would be observed close to the equatorial South Pole with respect to expectations for isotropy, while a decrease would be observed close to the edge of the directional exposure in the equatorial Northern Hemisphere. This would result in the observation of a fake dipole. A convenient way to exhibit this effect is to separate the dipole in two components: the component of the dipole in the equatorial plane r_⊥, and the component along Earth's rotation axis, r_∥. While r_⊥ is expected to be affected only by time-dependent effects, r_∥ is on the other hand the relevant quantity sensitive to time-independent effects such as the geomagnetic one.

Estimations of r_⊥ and r_∥ obtained by accounting or not for geomagnetic effects are given in Table 2, in two different energy ranges. These estimations are obtained from the recovered $\overline{a}_{1m}$ coefficients: $\overline{r}_\parallel =\sqrt{3}\overline{a}_{10}/\overline{a}_{00}$, and $\overline{r}_\perp =[3(\overline{a}_{11}^2+\overline{a}_{1-1}^2]^{0.5}/\overline{a}_{00}$. It can be seen that the main effect of the geomagnetic corrections is a shift in $\overline{r}_\parallel$ of about 1.2%. In the energy range 1 ⩽ E/EeV ⩽ 4, this shift is significant, $\overline{r}_\parallel$ changing from −2.2% to −1.0% with an uncertainty amounting to 0.4%. Above 4 EeV, the net correction is of the same order, though the statistical uncertainties are larger. In contrast, $\overline{r}_\perp$ remains unchanged in both cases, as expected.

Table 2. Influence of Shower Size Corrections for Geomagnetic Effects on the Component of the Dipole in the Equatorial Plane and on the One Along Earth's Rotation Axis

ΔE	$\overline{r}_\perp ^{{\rm uncorr}}$	$\overline{r}_\perp$	$\overline{r}_\parallel ^{{\rm uncorr}}$	$\overline{r}_\parallel$
(EeV)	(%)	(%)	(%)	(%)
1–4	0.9 ± 0.3	0.9 ± 0.3	−2.2 ± 0.4	−1.0 ± 0.4
>4	1.8 ± 1.0	2.1 ± 1.0	−4.1 ± 1.7	−3.0 ± 1.7

Download table as:  ASCIITypeset image

In this section, we study to which extent the procedure used to obtain the attenuation curve in Section 3.3 might influence the determination of the anisotropy parameters.

To convert the shower size into energy, we explained and applied in Section 3.3 the constant intensity cut method for showers with $S_{38^\circ }\ge 22$VEM, that is, just above the threshold energy for full efficiency. The value of the parameter a obtained in these conditions is consistent within the statistical uncertainties with the one previously reported when applying the same constant intensity cut method for showers with $S_{38^\circ }\ge 47$VEM. Opposite to this, the value obtained for the coefficient b differs by more than three standard deviations. Such a difference might be expected from both the evolution of the maximum of the showers and from an eventual change in composition with energy, but it may also be due to energy- and angle-dependent resolution effects mimicking a real evolution with energy.

With a different attenuation curve, some events would be reconstructed in the adjacent energy intervals to an extent which depends on the change of the attenuation curve with zenith angle. For that reason, the determination of anisotropy parameters might be altered by this effect.

Disentangling the real evolution of the attenuation curve from the energy from resolution effects is out of the scope of this paper and will be addressed elsewhere. Here, we restrict ourselves to probing the effect that a real energy dependence would have on the determination of anisotropy parameters. To do so, we choose to fit the values of the coefficient b obtained for $S_{38^\circ }=22$VEM and $S_{38^\circ }=47$VEM through a linear dependence with the logarithm of $S_{38^\circ }$. Below and above these values, the behavior of b(E) is obtained by extrapolating this energy dependence. In this way, the changes in the anisotropy parameters are probed in extreme conditions.

Repeating the whole chain of analysis with this new attenuation curve, it turns out that the reconstructed dipole parameters are only marginally affected by this change, as illustrated in the top and middle panels of Figure 10. Meanwhile, both reconstructed quadrupole amplitudes in the energy interval 2 ⩽ E/EeV ⩽ 4 are reduced in such a way that they lie now just below the 99% upper bounds for isotropy. Conversely, the amplitudes in the energy interval 1 ⩽ E/EeV ⩽ 2 are slightly increased. Below 4 EeV, the determination of the attenuation curve thus appears to bring some systematic uncertainties for determining the quadrupole amplitudes. The two extreme extrapolations performed in this analysis (i.e., b constant with the energy or linearly dependent with the logarithm of the energy) allow us to bracket the possible values.

Zoom In Zoom Out Reset image size

In Section 3, we presented the procedure adopted to account for the changes in shower size due to weather and geomagnetic effects. Since the coefficients α_P, α_ρ, and β_ρ in Equation (1) were extracted from real data, they suffer from statistical uncertainties which may impact in a systematic way the corrections made on S(1000), and consequently may also impact the anisotropy parameters derived from the data set. Besides, the determination of g₁ and g₂ in Equation (2) is based on the simulation of showers. Both the systematic uncertainties associated with the different interaction models and primary masses and the statistical uncertainties related to the procedure used to extract g₁ and g₂ constitute a source of systematic uncertainties on the anisotropy parameters.

To quantify these systematic uncertainties, we repeated the whole chain of analysis on a large number of modified data sets. Each modified data set is built by randomly sampling the coefficients α_P, α_ρ, and β_ρ (or g₁ and g₂ when dealing with geomagnetic effects) according to the corresponding uncertainties and correlations between parameters through the use of a Gaussian probability distribution function. For each new set of correction coefficients, new sets of anisotropy parameters are then obtained. The rms of each resulting distribution for each anisotropy parameter is the systematic uncertainty that we assign. Results are shown in Figure 10, in terms of the dipole and quadrupole amplitudes as a function of the energy. Balanced against the statistical uncertainties in the original analysis (shown by the bands), it is apparent that both sources of systematic uncertainties have a negligible impact on each reconstructed anisotropy amplitude.