Observation of a Large-scale Anisotropy in the Arrival Directions of Cosmic Rays above 8×1018 eV

Cosmic rays are atomic nuclei arriving from outer space that reach the highest energies observed in nature. Clues to their origin come from studying the distribution of their arrival directions. Using 3×104 cosmic rays above 8×1018 electron volts, recorded with the Pierre Auger Observatory from a total exposure of 76,800 square kilometers steradian year, we report an anisotropy in the arrival directions. The anisotropy, detected at more than the 5.2σ level of significance, can be described by a dipole with an amplitude of 6.5+1.3 −0.9% towards right ascension αd = 100± 10 degrees and declination δd = −24+12 −13 degrees. That direction indicates an extragalactic origin for these ultra-high energy particles. Particles with energies ranging from below 109 eV up to beyond 1020 eV, known as cosmic rays, constantly hit the Earth’s atmosphere. The flux of these particles steeply decreases as their energy increases; for energies above 10 EeV (1 EeV ≡ 1018 eV), the flux is about one particle per km2 per year. The existence of cosmic rays with such ultra-high energies has been known for more than 50 years [1, 2], but the sites and mechanisms of their production remain a mystery. Information about their origin can be obtained from the study of the energy spectrum and the mass composition of cosmic rays. However, the most direct evidence of the location of the progenitors is expected to come from studies of the distribution of their arrival directions. Indications of possible hot spots in arrival directions for cosmic rays with energy above 50 EeV have been reported by the Pierre Auger and Telescope Array Collaborations [3, 4], but the statistical significance of these results is low. We report the observation, significant at a level of more than 5.2σ, of a large-scale anisotropy in arrival directions of cosmic rays above 8 EeV. Above 1014 eV, cosmic rays entering the atmosphere create cascades of particles (called extensive air-showers) that are sufficiently large to reach the ground. At 10 EeV, an extensive airshower (hereafter shower) contains ∼1010 particles spread over an area of ∼20 km2 in a thin disc moving close to the speed of light. The showers contain an electromagnetic component (electrons, positrons and photons) and a muonic component that can be sampled using arrays of particle detectors. Charged particles in the shower also excite nitrogen molecules in the air, producing fluorescence light that can be observed with telescopes during clear nights. The Pierre Auger Observatory, located near the city of Malargüe, Argentina, at latitude 35.2◦S, is designed to detect showers produced by primary cosmic rays above 0.1 EeV. It is a hybrid system, a combination of an array of particle detectors and a set of telescopes used to detect the fluorescence light. Our analysis is based on data gathered from 1600 water-Cherenkov detectors deployed over an area of 3000 km2 on a hexagonal grid with 1500-m spacing. Each detector contains 12 tonnes of ultrapure water in a cylindrical container, 1.2 m deep and 10 m2 in area, viewed by three 9-inch photomultipliers. A full description of the observatory, together with details of the methods used to reconstruct the arrival directions and energies of events, has been published [5]. *correspondence to: auger spokespersons@fnal.gov †The authors with their affiliations appear at the end of this article. 1 ar X iv :1 70 9. 07 32 1v 1 [ as tr oph .H E ] 2 1 Se p 20 17 It is difficult to locate the sources of cosmic rays, as they are charged particles and thus interact with the magnetic fields in our galaxy and the intergalactic medium that lies between the sources and Earth. They undergo angular deflections with amplitude proportional to their atomic number Z, to the integral along the trajectory of the magnetic field (orthogonal to the direction of propagation), and to the inverse of their energy E. At E ≈ 10 EeV, the best estimates for the mass of the particles [6] lead to a mean value for Z between 1.7 and 5. The exact number derived is dependent on extrapolations of hadronic physics, which are poorly understood because they lie well beyond the observations made at the Large Hadron Collider. Magnetic fields are not well constrained by data, but if we adopt recent models of the Galactic magnetic field [7,8], typical values of the deflections of particles crossing the Galaxy are a few tens of degrees for E/Z = 10 EeV, depending on the direction considered [9]. Extragalactic magnetic fields may also be relevant for cosmic rays propagating through intergalactic space [10]. However, even if particles from individual sources are strongly deflected, it remains possible that anisotropies in the distribution of their arrival directions will be detectable on large angular scales, provided the sources have a nonuniform spatial distribution or, in the case of a single dominant source, if the cosmic-ray propagation is diffusive [11–14]. Searches for large-scale anisotropies are conventionally made by looking for nonuniformities in the distribution of events in right ascension [15,16] because, for arrays of detectors that operate with close to 100% efficiency, the total exposure as a function of this angle is almost constant. The nonuniformity of the detected cosmic-ray flux in declination (fig. S1) imprints a characteristic nonuniformity in the distribution of azimuth angles in the local coordinate system of the array. From this distribution it becomes possible to obtain information on the three components of a dipolar model. Event observations, selection, and calibration We analyzed data recorded at the Pierre Auger Observatory between 1 January 2004 and 31 August 2016, from a total exposure of about 76,800 km2 sr year. The 1.2-m depth of the waterCherenkov detectors enabled us to record events at a useful rate out to large values of the zenith angle, θ. We selected events with θ < 80◦ enabling the declination range −90◦ < δ < 45◦ to be explored, thus covering 85% of the sky. We adopted 4 EeV as the threshold for selection; above that energy, showers falling anywhere on the array are detected with 100% efficiency [17]. The arrival directions of cosmic rays were determined from the relative arrival times of the shower front at each of the triggered detectors; the angular resolution was better than 1◦ at the energies considered here [5]. Two methods of reconstruction have been used for showers with zenith angles above and below 60◦ [17, 18]. These have to account for the effects of the geomagnetic field [17, 19] and, in the case of showers with θ < 60◦, also for atmospheric effects [20] because systematic modulations to the rates could otherwise be induced (see supplementary materials). The energy estimators for both data sets were calibrated using events detected simultaneously by the water-Cherenkov detectors and the fluorescence telescopes, with a quasi-calorimetric determination of the energy coming from the fluorescence measurements. The statistical uncertainty in the energy determination is 16% above 4 EeV and 12% above 10 EeV, whereas the systematic uncertainty on the absolute energy scale, common to both data sets, is 14% [21]. Evidence that the analyses of the events with θ < 60◦ and of those with 60◦ < θ < 80◦ are consistent with each other comes from the energy spectra determined for the two angular bands. The spectra agree within the statistical uncertainties over the energy range of interest [22]. We consider events in two energy ranges, 4 EeV < E < 8 EeV and E ≥ 8 EeV, as adopted in previous analyses (e.g., [23–25]). The bin limits follow those chosen previously in [26, 27]. The median energies for these bins are 5.0 EeV and 11.5 EeV, respectively. In earlier work [23–25], the event selection required that the station with the highest signal be surrounded by six operational detectors—a demanding condition. The number of triggered stations is greater than four for 99.2% of all events above 4 EeV and for 99.9% of events above 8 EeV, making it possible to use

Particles with energies ranging from below 10 9 eV up to beyond 10 20 eV, known as cosmic rays, constantly hit the Earth's atmosphere. The flux of these particles steeply decreases as their energy increases; for energies above 10 EeV (1 EeV ≡ 10 18 eV), the flux is about one particle per km 2 per year. The existence of cosmic rays with such ultra-high energies has been known for more than 50 years [1,2], but the sites and mechanisms of their production remain a mystery. Information about their origin can be obtained from the study of the energy spectrum and the mass composition of cosmic rays. However, the most direct evidence of the location of the progenitors is expected to come from studies of the distribution of their arrival directions. Indications of possible hot spots in arrival directions for cosmic rays with energy above 50 EeV have been reported by the Pierre Auger and Telescope Array Collaborations [3,4], but the statistical significance of these results is low. We report the observation, significant at a level of more than 5.2σ, of a large-scale anisotropy in arrival directions of cosmic rays above 8 EeV.
Above 10 14 eV, cosmic rays entering the atmosphere create cascades of particles (called extensive air-showers) that are sufficiently large to reach the ground. At 10 EeV, an extensive airshower (hereafter shower) contains ∼10 10 particles spread over an area of ∼20 km 2 in a thin disc moving close to the speed of light. The showers contain an electromagnetic component (electrons, positrons and photons) and a muonic component that can be sampled using arrays of particle detectors. Charged particles in the shower also excite nitrogen molecules in the air, producing fluorescence light that can be observed with telescopes during clear nights.
The Pierre Auger Observatory, located near the city of Malargüe, Argentina, at latitude 35.2 • S, is designed to detect showers produced by primary cosmic rays above 0.1 EeV. It is a hybrid system, a combination of an array of particle detectors and a set of telescopes used to detect the fluorescence light. Our analysis is based on data gathered from 1600 water-Cherenkov detectors deployed over an area of 3000 km 2 on a hexagonal grid with 1500-m spacing. Each detector contains 12 tonnes of ultrapure water in a cylindrical container, 1.2 m deep and 10 m 2 in area, viewed by three 9-inch photomultipliers. A full description of the observatory, together with details of the methods used to reconstruct the arrival directions and energies of events, has been published [5].
It is difficult to locate the sources of cosmic rays, as they are charged particles and thus interact with the magnetic fields in our galaxy and the intergalactic medium that lies between the sources and Earth. They undergo angular deflections with amplitude proportional to their atomic number Z, to the integral along the trajectory of the magnetic field (orthogonal to the direction of propagation), and to the inverse of their energy E. At E ≈ 10 EeV, the best estimates for the mass of the particles [6] lead to a mean value for Z between 1.7 and 5. The exact number derived is dependent on extrapolations of hadronic physics, which are poorly understood because they lie well beyond the observations made at the Large Hadron Collider. Magnetic fields are not well constrained by data, but if we adopt recent models of the Galactic magnetic field [7,8], typical values of the deflections of particles crossing the Galaxy are a few tens of degrees for E/Z = 10 EeV, depending on the direction considered [9]. Extragalactic magnetic fields may also be relevant for cosmic rays propagating through intergalactic space [10]. However, even if particles from individual sources are strongly deflected, it remains possible that anisotropies in the distribution of their arrival directions will be detectable on large angular scales, provided the sources have a nonuniform spatial distribution or, in the case of a single dominant source, if the cosmic-ray propagation is diffusive [11][12][13][14].
Searches for large-scale anisotropies are conventionally made by looking for nonuniformities in the distribution of events in right ascension [15,16] because, for arrays of detectors that operate with close to 100% efficiency, the total exposure as a function of this angle is almost constant. The nonuniformity of the detected cosmic-ray flux in declination (fig. S1) imprints a characteristic nonuniformity in the distribution of azimuth angles in the local coordinate system of the array. From this distribution it becomes possible to obtain information on the three components of a dipolar model.

Event observations, selection, and calibration
We analyzed data recorded at the Pierre Auger Observatory between 1 January 2004 and 31 August 2016, from a total exposure of about 76,800 km 2 sr year. The 1.2-m depth of the water-Cherenkov detectors enabled us to record events at a useful rate out to large values of the zenith angle, θ. We selected events with θ < 80 • enabling the declination range −90 • < δ < 45 • to be explored, thus covering 85% of the sky. We adopted 4 EeV as the threshold for selection; above that energy, showers falling anywhere on the array are detected with 100% efficiency [17]. The arrival directions of cosmic rays were determined from the relative arrival times of the shower front at each of the triggered detectors; the angular resolution was better than 1 • at the energies considered here [5].
Two methods of reconstruction have been used for showers with zenith angles above and below 60 • [17,18]. These have to account for the effects of the geomagnetic field [17,19] and, in the case of showers with θ < 60 • , also for atmospheric effects [20] because systematic modulations to the rates could otherwise be induced (see supplementary materials). The energy estimators for both data sets were calibrated using events detected simultaneously by the water-Cherenkov detectors and the fluorescence telescopes, with a quasi-calorimetric determination of the energy coming from the fluorescence measurements. The statistical uncertainty in the energy determination is 16% above 4 EeV and 12% above 10 EeV, whereas the systematic uncertainty on the absolute energy scale, common to both data sets, is 14% [21]. Evidence that the analyses of the events with θ < 60 • and of those with 60 • < θ < 80 • are consistent with each other comes from the energy spectra determined for the two angular bands. The spectra agree within the statistical uncertainties over the energy range of interest [22].
We consider events in two energy ranges, 4 EeV < E < 8 EeV and E ≥ 8 EeV, as adopted in previous analyses (e.g., [23][24][25]). The bin limits follow those chosen previously in [26,27]. The median energies for these bins are 5.0 EeV and 11.5 EeV, respectively. In earlier work [23][24][25], the event selection required that the station with the highest signal be surrounded by six operational detectors-a demanding condition. The number of triggered stations is greater than four for 99.2% of all events above 4 EeV and for 99.9% of events above 8 EeV, making it possible to use events with only five active detectors around the one with the largest signal. With this more relaxed condition, the effective exposure is increased by 18.5%, and the total number of events increases correspondingly from 95,917 to 113,888. The reconstruction accuracy for the additional events is sufficient for our analysis (see supplementary materials and fig. S4).

Rayleigh analysis in right ascension
A standard approach for studying the large-scale anisotropies in the arrival directions of cosmic rays is to perform a harmonic analysis in right ascension, α. The first-harmonic Fourier components are given by The sums run over all N detected events, each with right ascension α i , with the normalization factor N = ∑ N i=1 w i . The weights, w i , are introduced to account for small nonuniformities in the exposure of the array in right ascension and for the effects of a tilt of the array towards the southeast (see supplementary materials). The average tilt between vertical and the normal to the plane on which the detectors are deployed is 0.2 • , so that the effective area of the array is slightly larger for showers arriving from the downhill direction. This introduces a harmonic dependence in azimuth of amplitude 0.3% × tan θ to the exposure. The effective aperture of the array is determined every minute. Because the exposure has been accumulated over more than 12 years, the total aperture is modulated by less than ∼0.6% as the zenith of the observatory moves in right ascension. Events are weighted by the inverse of the relative exposure to correct these effects ( fig. S2).
The amplitude r α and phase ϕ α of the first harmonic of the modulation are obtained from (2) Table 1 shows the harmonic amplitudes and phases for both energy ranges. The statistical uncertainties in the Fourier amplitudes are √ 2/N ; the uncertainties in the amplitude and phase correspond to the 68% confidence level of the marginalized probability distribution functions. The rightmost column shows the probabilities that amplitudes larger than those observed could arise by chance from fluctuations in an isotropic distribution. These probabilities are calculated as P(r α ) = exp(−N r 2 α /4) [28]. For the lower energy bin (4 EeV < E < 8 EeV), the result is consistent with isotropy, with a bound on the harmonic amplitude of <1.2% at the 95% confidence level. For the events with E ≥ 8 EeV, the amplitude of the first harmonic is 4.7 +0.8 −0.7 %, which has a probability of arising by chance of 2.6×10 −8 , equivalent to a two-sided Gaussian significance of 5.6σ. The evolution of the significance of this signal with time is shown in fig. S3; the dipole became more significant as the exposure increased. Allowing for a penalization factor of 2 to account for the fact that two energy bins were explored, the significance is reduced to 5.4σ. Further penalization for the four additional lower energy bins examined in [23] has a similarly mild impact on the significance, which falls to 5.2σ. The maximum of the modulation is at right ascension of 100 • ± 10 • . The maximum of the modulation for the 4 EeV < E < 8 EeV bin, at 80 • ± 60 • , is compatible with the one determined in the higher-energy bin, although it has high uncertainty and the amplitude is not statistically significant. Table S1 shows that results obtained under the stricter trigger condition and for the additional events gained after relaxing the trigger are entirely consistent with each other. Figure 1 shows the distribution of the normalized rate of events above 8 EeV as a function of right ascension. The sinusoidal function corresponds to the first harmonic; the distribution is compatible with a dipolar modulation: χ 2 /n = 10.5/10 for the first-harmonic curve and χ 2 /n = 45/12 for a constant function (where n is the number of degrees of freedom, equal to the number of points in the plot minus the number of parameters of the fit).   Table 1, which displays good agreement with the data (χ 2 /n = 10.5/10); the dashed line shows a constant function.
The distribution of events in equatorial coordinates, smoothed with a 45 • radius top-hat function to better display the large-scale features, is shown in Fig. 2.

Reconstruction of the three-dimensional dipole
In the presence of a three-dimensional dipole, the Rayleigh analysis in right ascension is sensitive only to its component orthogonal to the rotation axis of Earth, d ⊥ . A dipole component in the direction of the rotation axis of Earth, d z , induces no modulation of the flux in right ascension, but does so in the azimuthal distribution of the directions of arrival at the array. A non-vanishing value of d z leads to a sinusoidal modulation in azimuth with a maximum toward the northern or the southern direction.
To recover the three-dimensional dipole, we combine the first-harmonic analysis in right ascension with a similar one in the azimuthal angle ϕ, measured counterclockwise from the east. The relevant component, b ϕ , is given by an expression analogous to that in Eq. 1, but in terms of the azimuth of the arrival direction of the shower rather than in terms of the right ascension. The results are b ϕ = −0.013 ± 0.005 in the 4 EeV < E < 8 EeV bin and b ϕ = −0.014 ± 0.008 in the E ≥ 8 EeV bin. The probabilities that larger or equal absolute values for b ϕ arise from an isotropic distribution are 0.8% and 8%, respectively.
Under the assumption that the dominant cosmic-ray anisotropy is dipolar, based on previous  100 ± 10 studies that found that the effects of higher-order multipoles are not significant in this energy range [25,29,30], the dipole components and its direction in equatorial coordinates (α d , δ d ) can be estimated from [25], where cos δ is the mean cosine of the declinations of the events, sin θ is the mean sine of the zenith angles of the events, and obs −35.2 • is the average latitude of the Observatory. For our data set, we find cos δ = 0.78 and sin θ = 0.65. The parameters describing the direction of the three-dimensional dipole are summarized in Table 2. For 4 EeV < E < 8 EeV, the dipole amplitude is d = 2.5 +1.0 −0.7 %, pointing close to the celestial south pole, at (α d , δ d ) = (80 • , −75 • ), although the amplitude is not statistically significant. For energies above 8 EeV, the total dipole amplitude is d = 6.5 +1.3 −0.9 %, pointing toward (α d , δ d ) = (100 • , −24 • ). In Galactic coordinates, the direction of this dipole is ( , b) = (233 • , −13 • ). This dipolar pattern is clearly seen in the flux map in Fig. 2. To establish whether the departures from a perfect dipole are just statistical fluctuations or indicate the presence of additional structures at smaller angular scales would require at least twice as many events.

Implications for the origin of high-energy cosmic rays
The anisotropy we have found should be seen in the context of related results at lower energies. Above a few PeV, the steepening of the cosmic ray energy spectrum has been interpreted as being due to efficient escape of particles from the Galaxy and/or because of the inability of the sources to accelerate cosmic rays beyond a maximum value of E/Z. The origin of the particles remains unknown. Although supernova remnants are often discussed as sources, evidence has been reported for a source in the Galactic center capable of accelerating particles to PeV energies [31]. Diffusive escape from the Galaxy is expected to lead to a dipolar component with a maximum near the Galactic center direction [32]. This is compatible with results obtained in the 10 15 to 10 18 eV range [15,16,23,24,33], which provide values for the phase in right ascension close to that of the Galactic center, α GC = 266 • .
Models proposing a Galactic origin up to the highest observed energies [34,35] are in increasing tension with observations. If the Galactic sources postulated to accelerate cosmic rays above EeV energies, such as short gamma-ray bursts or hypernovae, were distributed in the disk of the Galaxy, a dipolar component of anisotropy is predicted with an amplitude that exceeds existing bounds at EeV energies [24,33]. In this sense, the constraint obtained here on the dipole amplitude (Table 2) for 4 EeV < E < 8 EeV further disfavors a predominantly Galactic origin. This tension could be alleviated if cosmic rays at a few EeV were dominated by heavy nuclei such as iron, but this would be in disagreement with the lighter composition inferred observationally at these energies [6]. The maximum of the flux might be expected to lie close to the Galactic center region, whereas the direction of the three-dimensional dipole determined above 8 EeV lies ∼125 • from the Galactic center. This suggests that the anisotropy observed above 8 EeV is better explained in terms of an extragalactic origin. Above 40 EeV, where the propagation should become less diffusive, there are no indications of anisotropies associated with either the Galactic center or the Galactic plane [36].
There have been many efforts to interpret the properties of ultrahigh-energy cosmic rays in terms of extragalactic sources. Because of Liouville's theorem, the distribution of cosmic rays must be anisotropic outside of the Galaxy for an anisotropy to be observed at Earth. An anisotropy cannot arise through deflections of an originally isotropic flux by a magnetic field. One prediction of anisotropy comes from the Compton-Getting effect [37], which results from the proper motion of the Earth in the rest frame of cosmic-ray sources, but the amplitude is expected to be only 0.6% [38], well below what has been observed. Other studies have predicted larger anisotropies. These assume that ultrahigh-energy cosmic rays originate from an inhomogeneous distribution of sources [13,14,39], or that they arise from a dominant source and then diffuse through intergalactic magnetic fields [11][12][13][14]. The resulting dipole amplitudes are predicted to grow with energy, reaching 5 to 20% at 10 EeV. These amplitudes depend on the cosmicray composition as well as the details of the source distribution. On average, the predictions are smaller for larger source densities or for more isotropically distributed sources. If the sources were distributed like galaxies, the distribution of which has a significant dipolar component [40], a dipolar cosmic-ray anisotropy would be expected in a direction similar to that of the dipole associated with the galaxies. This effect would be due to the excess of cosmic-ray sources in this direction and is different from the Compton-Getting effect due to the Earth's motion with respect to the rest frame of cosmic rays. For the infrared-detected galaxies in the 2MRS catalogue [40], the flux-weighted dipole points in Galactic coordinates in the direction ( , b = (251 • , 38 • ). In this coordinate system, the dipole we detect for cosmic rays above 8 EeV is in the direction (233 • , −13 • ), about 55 • away from that of the 2MRS dipole.
For an extragalactic origin, the Galactic magnetic fields modify the direction of the dipole observed at Earth relative to its direction outside the Galaxy. For illustration, Fig. 3 shows a map of the flux above 8 EeV in which the direction of the cosmic-ray dipole is shown along with the direction towards the 2MRS dipole. The arrows in the plot indicate how a dipolar distribution of cosmic rays, in the same direction as the 2MRS dipole outside the Galaxy, would be affected by the Galactic magnetic field [8]. average values for Z ∼ 1.7 to 5 at 10 EeV, these represent typical values of E/Z for the cosmic rays contributing to the observed dipole. The agreement between the directions of the dipoles is improved by adopting these assumptions about the charge composition and the deflections in the Galactic magnetic field. For these directions, the deflections within the Galaxy will also lead to a lowering of the amplitude of the dipole to about 90% and 70% of the original value, for E/Z = 5 EeV and 2 EeV, respectively. The lower amplitude in the 4 EeV < E < 8 EeV bin might also be the result of stronger magnetic deflections at lower energies. Our findings constitute the observation of an anisotropy in the arrival direction of cosmic rays with energies above 8 EeV. The anisotropy can be well represented by a dipole with an amplitude of 6.5 +1.3 −0.9 % in the direction of right ascension α d = 100 ± 10 • and declination δ d = −24 +12 −13 • . By comparing our results with phenomenological predictions, we find that the magnitude and direction of the anisotropy support the hypothesis of an extragalactic origin for the highest-energy cosmic rays, rather than sources within the Galaxy.

Acknowledgments
The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment from the technical and administrative staff in Malargüe, and the financial support from a number of funding agencies in the participating countries. The full acknowledgments are in the Supplementary Materials. The Pierre Auger Collaboration will make public the data reproduced in the Figures of this paper on the Auger website on the link to this publication: www.auger.org/data/science2017.tar.gz

Event weights
The Fourier harmonic analysis we perform accounts for modulations in the exposure (of instrumental origin) as well as the effects due to the tilt of the surface detector array. To achieve this each event is weighted by a factor where θ i and φ i are the zenith and azimuth angles of the events and α 0 i is the right ascension directly above the array at the time the i-th event was detected. The term enclosed by the round brackets accounts for the effects of the slight tilt of the array. The average inclination from the vertical is θ tilt ≈ 0.2 • , in the direction φ tilt ≈ −30 • , i.e. 30 • South of the Easterly direction. The tilt term affects only the Fourier analysis in azimuth, and thus the dipole component d z . The second term in Eq. (S1), ∆N cell (α 0 ), allows for the fact that the effective aperture of the observatory is not uniform in sidereal time. This factor corresponds to the relative number of detector cells, i.e. the active detectors surrounded by at least five other active detectors, present when the right ascension of the zenith of the observatory equals α 0 within binning accuracy. It is obtained by adding the number of cells over the whole period of observations, with dead times due to power failures or to communication or acquisition problems discarded. The total number of cells within each α 0 bin is normalized to the average value [22]. ∆N cell (α 0 ) is plotted in Fig. S2 with a bin width of 1 • . This term affects only the Fourier analysis in right ascension, and thus the dipole component d ⊥ . After more than 12 years of continuous operation of the observatory the normalized number of cells shows variations that are smaller than ±0.6%. If the effects of the modulation in the number of cells were not taken into account through the weights, a spurious contribution to d ⊥ of amplitude 0.05% in the direction α ≈ 145 • would be induced. This contribution is an order of magnitude smaller than the statistical uncertainty in the determination of this dipole component, having then a marginal effect. If the effects of the tilt were not taken into account, a spurious contribution to d z of −0.4% would be induced. Had we not introduced weighting to the Fourier analysis we would have obtained: for 4 EeV < E < 8 EeV, a total dipole amplitude d = 2.

Energy reconstruction
Two methods of reconstruction have been used for showers with zenith angles above and below 60 • [16,17]. The energy estimator used for showers with θ < 60 • is the signal reconstructed at 1000 m from the shower core. This signal is corrected for atmospheric effects [18] that would otherwise introduce systematic modulations to the rates as a function of time of day or season. This could result in spurious influences on the distribution in sidereal time (a time scale that is based on the Earth's rate of rotation measured relative to the fixed stars rather than the Sun, corresponding to 366.25 cycles/year) and hence could be a source of systematic effects for the anisotropies inferred. The atmospheric effects arise from the dependences of the longitudinal and lateral attenuation of the electromagnetic component of air showers on atmospheric conditions, in particular temperature and pressure. If not corrected, these could cause a modulation of the rates of up to ±1.7% in solar time. The energy estimator is also corrected for geomagnetic effects [19] as otherwise a systematic modulation of amplitude ∼0.7% would be induced in the azimuthal distributions.
The particles arriving at the ground in showers with θ > 60 • are predominantly muons. As the atmospheric thickness traversed by a shower is proportional to sec θ, at those zenith angles θ < 60 60 < θ < 80 θ < 80 Figure S1: Exposure as a function of declination. The separate contributions from events with θ < 60 • and 60 • < θ < 80 • are also displayed.
the electromagnetic component is almost completely absorbed so that atmospheric effects are negligible. For these large angles the energy estimator is based on the muon content relative to that in a simulated proton-shower of 10 EeV, with the geomagnetic deflections of muons accounted for in the reconstruction of these events [17]. After applying the corrections to those events with θ < 60 • , the amplitude of modulation in solar time (365.25 cycles/year) for the whole data set (with θ < 80 • and E > 4 EeV) is reduced to 0.5 ± 0.4%. This is obtained from the first harmonic Fourier analysis of the arrival times as a function of the hour of the day. The residual effect in right ascension, after averaging over more than 12 years, is less than one part in a thousand. As a further check, the amplitude of modulation at the anti-sidereal frequency (364.25 cycles/year) is 0.5 ± 0.4%, consistent with the absence of residual systematic effects. The results of the solar and anti-sidereal amplitudes in the two separate energy bins are also consistent with the absence of systematic effects, being, for 4 EeV < E < 8 EeV of 0.6 ± 0.5% in solar and of 0.4 ± 0.5% in anti-sidereal, while for E > 8 EeV they are 0.7 ± 0.8% in solar and of 1.1 ± 0.8% in anti-sidereal.

Accuracy of the reconstruction of events obtained with relaxed trigger
Taking events passing the stricter cuts (with six active detectors surrounding the one with the highest signal) and re-analyzing them after removing one of the six detectors, we find that for E ≥ 8 EeV (4 EeV < E < 8 EeV) the difference between the reconstructed directions has an average of 0.3 • (0.4 • ), with 90% of the events having an angular difference smaller than 0.7 • (1.2 • ). The energy estimates differ on average by only 0.2% (0.3%), with a dispersion of ∼5% (8%). These differences are well below the experimental uncertainties of these two parameters. The distribution of the differences in arrival directions and reconstructed energies between the original event satisfying the strict trigger and those with a missing detector around the one with the highest signal are shown in fig. S4. Events passing the stricter cut from two years of data (2013 and 2014) were analyzed, leading to a total of artificial events passing the relaxed trigger of 65,000 for 4 EeV < E < 8 EeV and 27,000 for E ≥ 8,EeV. The best-fitting first and second harmonics are overlaid. The first harmonic has an amplitude of (0.06 ± 0.02)% and the second harmonic has an amplitude of (0.15 ± 0.02)%.

Tight Triggers Relaxed Triggers
All Triggers Figure S3: Evolution of the probability that the signal arises by chance as a function of time. The evolution of the probability that in an isotropic distribution the first-harmonic amplitude in right ascension be larger or equal than the one measured is shown for events with E ≥ 8 EeV. This is plotted as a function of the exposure accumulated over the years. The solid blue line corresponds to the signal from the combination of the two triggers while the black and red lines refer to data from the tight and relaxed triggers respectively (Table S1). The values corresponding to 5σ and 5.5σ are indicated as horizontal lines.  Figure S4: The difference between measurements of energy and direction for the strict and relaxed trigger conditions. Top panel: Distribution of relative difference between the reconstructed energies of the event satisfying the strict trigger condition and those obtained by artificially removing one of the six detectors surrounding the one with the largest signal. Bottom panel: corresponding distribution of the angular separation between the reconstructed arrival directions. For E ≥ 8 EeV (4 EeV < E < 8 EeV) the difference between the reconstructed directions has an average of 0.3 • (0.4 • ), with 90% of the events having an angular difference smaller than 0.7 • (1.2 • ). The energy estimates differ on average by only 0.2% (0.3%), with a dispersion of ∼5% (8%).