Introduction

In the past decade, significant progress has been made in the search for the origin of ultra-high-energy cosmic rays (UHECRs) thanks to the data from the Pierre Auger Observatory (Pierre Auger Collaboration 2015b) and Telescope Array (Telescope Array Collaboration 2013). The study of the arrival directions is a key element in gaining insight into the possible sources of UHECRs. At small and intermediate angular scales, interesting hints of anisotropies have been reported, the most significant being the overdensity in the Centaurus region, with a significance of 4σ (see Pierre Auger Collaboration 2022, 2023d for the latest results). Regarding large angular scales, a dipolar modulation in R.A. at energies above 8 EeV has been determined with a significance above 5σ (Pierre Auger Collaboration 2017c). The direction of this dipole, lying ∼115°away from the Galactic center, suggests an extragalactic origin for the cosmic rays above this energy threshold. 97 Furthermore, an approximately linear growth of the dipole amplitude as a function of energy has been measured (Pierre Auger Collaboration 2018). This could be a consequence of the energy losses that cosmic rays suffer in interactions with the background radiation, which lead to a decrease of the horizon of cosmic-ray sources as the energy increases (K. Greisen 1966; G. T. Zatsepin & V. A. Kuz'min 1966). This shrinking of the horizon is expected to induce an increase in the dipole amplitude as a function of energy, due to the growing relative contribution of nearby sources, whose distribution is more inhomogeneous. In the alternative scenario in which only one or a few nearby sources give the dominant contribution to the flux above 4 EeV, a rising dipolar amplitude with energy is also expected owing to the growth of the UHECR magnetic rigidity, which implies smaller deflections of the particles in the intervening magnetic fields. 98 For energies below 4 EeV, where only the anisotropies in R. A. can be studied, small amplitudes of the equatorial dipole, below the 1% level, have been measured, which are compatible with isotropic expectations for the present accumulated statistics (Pierre Auger Collaboration 2020a). However, a change in the R.A. phase at energies below 4 EeV is apparent, with the observed values being close to the R.A. of the Galactic center. This change could be indicative of a transition in the origin of the anisotropies from an extragalactic one at energies above a few EeV to a Galactic one at lower energies, or alternatively due to the effects of the Galactic magnetic field on an extragalactic flux (S. Mollerach et al. 2022).

In this work, to gain further insight into these results, the statistics of 19 yr of data from the Pierre Auger Observatory are analyzed, which represents an increase of ∼55% with respect to Pierre Auger Collaboration (2018) and of ∼30% with respect to Pierre Auger Collaboration (2020a).

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 97 Due to the smaller statistics of Telescope Array (∼10 times smaller), a large-scale dipolar anisotropy has not been confirmed in their data set, with a 99% CL upper limit on the first-harmonic amplitude of 7.3% (Telescope Array Collaboration 2020). 98 For the source model that best fits the spectrum and composition results obtained in Pierre Auger Collaboration (2023a), the cosmic-ray horizon is shrunk from a few hundreds of megaparsecs to a few tens of megaparsecs between 10 and 100 EeV, and the mean rigidity changes from ∼4 to ∼8 EV in the same energy range.

An update on the angular power spectrum (Pierre Auger Collaboration 2017a) is also presented to determine whether higher multipoles are present. Note that even a pure dipolar distribution at the boundaries of the galaxy could be transformed, due to the deflections of cosmic rays in the Galactic magnetic field, into a dipolar distribution plus higherorder multipoles.

Finally, predictions are included, based on our combined fit of spectrum and composition data (Pierre Auger Collaboration 2023a), for the dipolar amplitude and direction and for the quadrupolar amplitude.

The Pierre Auger Observatory and the Data Sets

The Pierre Auger Observatory (Pierre Auger Collaboration 2015b) is the largest cosmic-ray observatory in the world, covering an area of 3000 km 2 . It is located near the city of Malargüe, Argentina (35.°2 south, 69.°5 west, 1400 m above sea level). It is a hybrid detector, consisting of a surface detector (SD) with 1660 water-Cherenkov stations and a fluorescence detector (FD) with 27 telescopes overlooking the array. The duty cycle of the SD is ∼100%, while that corresponding to the FD, which operates on clear moonless nights, is ∼15%. The main advantage of having both detectors is that the quasi-calorimetric energy determination obtained with the FD can be used to calibrate the energy estimate of the events registered with the SD, using the hybrid events detected simultaneously by the SD and the FD.

The main surface array (SD-1500) consists of 1600 detectors on a triangular grid with a spacing of 1500 m. There is a denser "infilled array" (SD-750) of 60 detectors with a spacing of 750 m, used to register events down to an energy threshold of ∼0.03 EeV. There is also a smaller array with a spacing of 433 m, but the data set has smaller statistics than the SD-750 and is not included in this work. The events recorded with the SD-1500 that are considered here were registered between 2004 January 1 and 2022 December 31. For the transitional years of 2021 and 2022, when the AugerPrime (Pierre Auger Collaboration 2023b) installation was underway, only those detectors in which the electronics had not been updated are used (resulting in an equivalent ∼1.6 yr of exposure instead of 2 yr).

The analyses above 4 EeV are made considering events with zenith angle θ < 80°, achieving an 85% coverage of the sky. The events are selected if the SD station with the largest signal is surrounded by at least five active stations and if the reconstructed core of the shower falls within an isosceles triangle of nearby active stations. To accurately compute the exposure as a function of time, those events that were registered during periods of unreliable data acquisition are removed, resulting in a total exposure of 123,000 km 2 sr yr. There are actually different reconstructions for the events with θ < 60°("vertical") and 60°< θ < 80°("inclined"). Full efficiency is reached above 2.5 EeV for the former events and above 4 EeV for the latter ones.

For studies using the SD-1500 array between 0.25 and 4 EeV only vertical events are considered (leading to a 71% coverage of the sky) and a stricter quality cut, requiring that the SD station with the largest signal should be surrounded by six active stations, is used. This data set has an exposure of 81,000 km 2 sr yr.

Events with energies down to 0.03 EeV with the SD-750 array are also used. In this data set, not only events registered using the standard station-level trigger algorithms are considered but also, for the first time, those detected with two other triggers introduced in mid-2013 to enhance the sensitivity of the array to small signals (Pierre Auger Collaboration 2021). Thanks to these triggers, full efficiency for the SD-750 array is reached above 0.2 EeV for events with θ < 55°. The accumulated exposure of the SD-750 array with θ < 55°between 2014 January 1 and 2021 December 13 is 269 km 2 sr yr.

The statistical uncertainty of the energy is ∼7% for E > 10 EeV and can be up to ∼20% for E ∼ 0.1 EeV, while the systematic uncertainty of the absolute energy scale is 14% (Pierre Auger Collaboration 2020b, 2021). The events have an angular resolution better than 0°. 9 for E > 10 EeV, which can degrade to 2°at lower energies (Pierre Auger Collaboration 2015b, 2023c).

The energies of vertical events are corrected for atmospheric (Pierre Auger Collaboration 2017b) and geomagnetic effects (Pierre Auger Collaboration 2011), so as to avoid spurious modulations in the distributions in R.A. and azimuth, respectively. For inclined events, the air shower cascades near ground level are composed mostly of muons, and the atmospheric effects are negligible for them, while the geomagnetic effects are directly included in the reconstruction.

Description of the Large Angular Scale Analyses

In this section, the methods used in the present work are briefly discussed. For further details see Pierre Auger Collaboration (2015a, 2017a, 2018, 2020a).

3D Dipole above 4 EeV

To reconstruct the 3D dipole, a separate Fourier analyses in R.A. (α) and azimuth (f) is performed. The amplitude r k

x and phase k x j of the event rate modulation are given by

where the harmonic amplitudes of order k are

with x = α or f and k = 1 for the dipolar component. N is the total number of events, w i are the weight factors, and

is the normalization factor. The weight factors w i are computed as

where ΔN cell (α 0 (t i )) is the relative variation of the total number of active detector cells for a given R.A. of the zenith of the observatory α 0 evaluated at the time t i at which the ith event is detected, and θ i and f i are the zenith and azimuth of the event, respectively. These weights are used to take into account the slightly nonuniform exposure of the array over time (due to the deployment of the array and the sporadic downtime of the detectors) and the tilt of the array (which is on average tilted 0°. 2 toward f 0 = -30°, with f 0 measured counterclockwise from the east). Both effects could induce spurious modulations if they are not accounted for.

The probability that an amplitude equal to or larger than r k a arises from a fluctuation from an isotropic distribution of events is given by ( ) ( ( ) Linsley 1975). Assuming that the dominant anisotropy in the distribution of arrival directions of cosmic rays is given by the dipolar component, d, the flux distribution can then be written as

where û is the arrival direction of the cosmic rays. The equatorial amplitude of the dipole d ⊥ , the north-south component d z , and the dipole direction in equatorial coordinates (α d , δ d ) are related to the first-harmonic amplitudes in R. A. and azimuth through ( )

is the mean cosine of the decl. of the events, sin 0.6525 q á ñ » is the mean sine of the zenith of the events, and l obs = -35°.2 is the latitude of the observatory.

3D Dipole and Quadrupole above 4 EeV

In the case in which a quadrupolar distribution is also present, the flux anisotropy can be parameterized as

where Q ij are the components of the quadrupole tensor (five of them are independent). The details on how the dipole and quadrupole components are estimated from the k = 1 and 2 harmonic amplitudes are described in Pierre Auger Collaboration (2015a).

Angular Power Spectrum Above 4 EeV

To search for anisotropies across various angular scales, it is convenient to decompose the distribution of observed events per unit solid angle

where ( ) W u ^is the relative exposure of the observatory, N is the total number of observed events, and

. The angular power spectrum of the scalar field

, where the coefficients a ℓm , encoding any anisotropy signature present in data, are derived from the multipolar expansion of

. Due to the incomplete sky coverage of the Pierre Auger Observatory, the estimation of the individual a ℓm coefficients cannot be carried out with relevant resolution as soon as ℓ 2 max > , and the same is true for the power spectrum (full-sky analyses are carried out by the Pierre Auger and Telescope Array Collaborations together; Pierre Auger & Telescope Array Collaborations 2014, 2023). However, it is possible to reconstruct the angular power spectrum within a statistical resolution independent of the bound ℓ max (O. Deligny et al. 2004), if the observed distribution of arrival directions represents a particular realization of an underlying stochastic process in which the anisotropies cancel in the ensemble average ( ( ˆ) ) u 0 áD ñ = and the second-order moment ( ( ˆ) (ˆ) u u áD D ¢ ñ) depends only on the angular separation between û and û¢. For a detailed discussion of the implications and restrictions of this hypothesis, see comments in the References section (Pierre Auger Collaboration 2017a). Under this hypothesis, the ensemble-average expectations of the power spectrum C ℓ á ñ and the "pseudo"-power spectrum Cℓ á ñ are related through

with the "pseudo"-power spectrum ˜| ˜| (

. The operator M, describing the mixing between the modes, is entirely determined by the knowledge of the relative coverage function (O. Deligny et al. 2004), and the f N 4 1 2 p term corresponds to an irreducible noise induced by Poisson fluctuations. Thus, for a cosmic-ray data set with a "pseudopower" spectrum Cℓ , the angular power spectrum Ĉℓ ,foot_1 which is an unbiased estimator in the ensemble average, is computed by

Modulation in R.A. Down to 0.03 EeV

When the trigger efficiency of the array is small, the systematic effects on the azimuth distribution cannot be completely corrected for, in particular the dependence of the trigger efficiency on azimuth due to the geometry of the layout of the array. Thus, for lower energies we restrict the large-scale studies to the modulation in R.A., which is proportional to the equatorial dipole component as shown in Equation (5).

Due to Earth's rotation, the exposure of the observatory is almost uniform in R.A. and above full efficiency the sources of spurious modulation (sporadic detector downtime and atmospheric effects) can be corrected for as described above. The effect of the tilt of the array is not relevant for the Fourier analysis in R.A. (it only affects the distribution in azimuth); thus, it is not included in the weights for these analyses.

Trigger effects at energies where the efficiency is low are difficult to control down to the required precision for reconstructing per-mille-level anisotropies. However, a method suitably constructed to be insensitive to such effects, called the East-West method (R. Bonino et al. 2011), can be used. This method, which is based on the difference between the counting rates of the events detected from the east sector and those from the west sector, is less sensitive than the Fourier method, but the systematics are under better control. The Fourier coefficients for the East-West method are

where α 0 (t i ) is the R.A. of the zenith of the array, and ξ i = 0 for events coming from the east (-π/2 < f < π/2) and ξ i = π for those coming from the west (π/2 < f < 3π/2). The Fourier amplitude, r, and phase, j, are related to the East-West method values, r EW and j EW , through

as in R. Bonino et al. (2011). As in the Fourier method, the probability of obtaining an amplitude larger than that expected from an isotropic distribution is ( ) ( ( )

. For the 0.25-0.5 EeV energy bin, one can also use the data of the SD-750 array with the standard Fourier method since the SD-750 array is fully efficient at these energies. The disadvantage of having lower statistics with the SD-750 array compared to those of the SD-1500 array at these energies is compensated by the fact that the number of events needed to have a statistical uncertainty with the Fourier method equal to that with the East-West method is approximately a factor four smaller (the statistical uncertainty obtained with the East-West method is larger by a factor cos 2 sin 2.1 p d q á ñ á ñ ~for θ < 60°; R. Bonino et al. 2011). For the energy bins below 0.25 EeV, the statistics of the SD-750 array is larger than that of the main array, and the East-West method is applied given that for those energies the trigger of the SD-750 array is not fully efficient. For the analyses in R.A., the geomagnetic corrections are not necessary (they are only relevant for the analyses in azimuth); therefore, these corrections are not done for the SD-750 data set. However, the atmospheric corrections in the reconstructed energy and the not completely uniform exposure in R.A. of the array are accounted for in the analyses with the SD-750 data set, like in the SD-1500 data set.

Results

3D Dipole above 4 EeV

The results for the 3D dipole reconstruction above full efficiency are listed in Table 1. For E 8 EeV, the significance of the dipolar modulation in R.A. is now at 6.8σ, and its significance in the 8-16 EeV energy bin is 5.7σ. The uncertainties reported correspond to the 68% confidence level (CL) of the respective marginalized probability distribution function.

In Figure 1, the flux above 8 EeV in equatorial coordinates and the distribution in R.A. of the rates of events (normalized to unity) for that energy threshold are depicted. In Figure 2 sky maps in Galactic coordinates are included showing the ratio between the number of observed events and the number expected from an isotropic distribution of arrival directions, for the energy bins of (4-8, 8-16, 16-32, 32) EeV. In Figure 2(d), the overdensity of events in the Centaurus region, (l, b) = (310°, 19°), is visible, contributing in the region corresponding to the dipole excess.

In Figure 3, the evolution with energy of the dipole direction and amplitude is plotted. A fit to the amplitude as a function of (2018). This growth in amplitude is possibly due to the larger relative contribution from the nearby sources for increasing energies, whose distribution is more inhomogeneous. Another effect that also results in an increase of the dipole amplitude with energy, secondary to this one, is the growth of mean rigidity of the particles, leading to smaller deflections and thus larger dipolar amplitudes (see Section 5 for a comparison of model predictions to data).

To show that there are no significant unaccounted-for systematics effects, in Table 2 the solar and antisidereal firstharmonic amplitudes (r 1 solar and r 1 antis ) and the probability of getting a larger amplitude from fluctuations of an isotropic distribution ( ( ) P r 1 solar  and ( ) P r 1 antis



) are reported, for each energy bin in which the Fourier analysis in R.A. was done at the sidereal frequency. The results are compatible with no significant modulation being present at these frequencies. Furthermore, for the Fourier analysis in azimuth, the a 1 f parameters, which measure the east-west modulation, were verified to be compatible with zero (since any flux modulation from the sky should be averaged owing to Earth's rotation).

The statistics of the E 8 EeV cumulative energy bin is large enough that the data set can be separated into timeordered subsets and tested for their stability. The results, dividing the data set into two (the date corresponding to half the data set is 2014 August 19) and four subsets (the dates corresponding to the end of each quarter of the data set being 2010 September 15, 2014 August 19, 2018 September 5, and 2022 December 31), are listed in Table 3. It is seen that the equatorial dipole amplitude and phase are consistent within the statistical uncertainties in the different time periods considered. Given that no time variation of the equatorial dipole modulation is found, an upper limit to its long-term rate of change of 0.003 yr -1 at the 95% CL can be set using a linear fit to the split data, which illustrates the long-term stability of the detector. Another check that was performed was to separate the data set into two zenith angle ranges, larger or smaller than 60°, corresponding to the different reconstructions (vertical and inclined), and the results are compatible with each other.

3D Dipole and Quadrupole above 4 EeV

In Table 4, the results obtained if a quadrupolar component is included in the Fourier analysis for each energy bin are shown. The dipole components d = (d x , d y , d z ), the five independent quadrupolar components Q ij , the quadrupole amplitude

, and the 99% CL upper limits Q UL are listed. It is seen that the quadrupolar components are not significant and that the dipole components are consistent with the results we obtain assuming only a dipole.

Angular Power Spectrum above 4 EeV

The measured power spectra for different energy bins are presented in Figure 4. The gray band, obtained from Monte Carlo simulations, corresponds to the 99% CL interval that would result from fluctuations of isotropic distributions. This band was by performing 50,000 isotropic simulations, each containing the same number of events as the data. The Ĉ1 and Ĉ2 obtained for all energy bins are consistent with the results from the previous Sections 4.1 and 4.2 using a Fourier analysis. Besides the significant dipolar pattern corresponding to Ĉ1 , the only Ĉℓ values that stand above the 99% CL of isotropic fluctuations are Ĉ17 , corresponding to an angular scale of ∼180°/ℓ ≈ 11°, and Ĉ8 , corresponding to an angular scale of ∼23°, for the energy bins of (4, 8) EeV and (16, 32) EeV, respectively. After statistical penalization for searches over different multipoles and energy bins (four independent energy bins × 20 multipole measurements = 80), the probabilities of these results arising from fluctuations of isotropy are 3.3% and 26.5%, respectively. All other Ĉℓ values in different energy bins are not significant.

The red lines indicate the upper limits on multipole amplitudes with 99% CL. They were obtained by using the following approach. For each possible amplitude C ℓ , the probability of reconstructing Ĉℓ rec given a true value of C ℓ ,

rec is estimated from simulations. For this, cosmic-ray skies were simulated by drawing events according to an underlying anisotropic distribution corresponding to C ℓ , considering an observatory with relative exposure W. The underlying anisotropic distribution is achieved by generating a ℓm coefficients uniformly distributed on the surface of a (2ℓ + 1)-dimensional hypersphere of radius (

) ℓ C 2 1 ℓ + . The reconstructed Ĉℓ rec for each simulated cosmic-ray sky are obtained using Equation (9). The amplitude C ℓ UL such that dC p C ℓ rec ℓ,data ò ¥ ( ˆ) C C , C L ℓ ℓ rec UL = is the relevant upper limit. 100 4.4. Modulation in R.A. Down to 0.03 EeV

The modulation in R.A. is studied for energies below 4 EeV using the data of the SD-1500 and SD-750 arrays, with the East-West method or a Fourier analysis, as described in Section 3.4. The results of these analyses, the equatorial amplitude (d ⊥ ), phase (α d ), and isotropic probability Note. The first two rows correspond to separating the data set into two timeordered bins, while the following rows correspond to separating the data set into four time-ordered bins.

100 Since this procedure can lead to upper limits tighter than the upper bounds obtained from isotropy with 99% CL, Ĉℓ,iso 99 , when the measured values of Ĉℓ,data are smaller than the expected average for isotropy, the upper limits presented are defined as max

in order to cope with this undesired behavior.

( 5, where it is indicated which method was used in each energy bin. The 99% CL upper limit, d UL ^, is reported. These upper limits are derived from the distribution for a dipolar anisotropy of unknown amplitude, marginalized over the dipole phase, as described in Pierre Auger Collaboration (2020a). In Figure 5, the results are shown, also including those obtained at lower energies by the IceCube and KASCADE-Grande Collaborations (IceCube Collaboration 2012, 2016; KASCADE-Grande Collaboration 2019).

Due to the small amplitudes (<1%) and the current amount of statistics, the results listed in Table 5 are not statistically significant ( (

), and it is not yet possible to draw firm conclusions. Nonetheless, it is suggestive that the amplitudes of the equatorial dipole increase as a function of energy, from below 1% to above 10%, and that the phases shift from close to the Galactic center to the opposite direction. These changes in amplitude and phase could suggest a transition of the origin of the anisotropies from Galactic cosmic rays to extragalactic ones. Another explanation for these changes could be the impact of the Galactic magnetic field on an extragalactic flux, as discussed in Section 5.

Model Predictions and Discussion

Possible explanations of the origin of the dipole measured above 8 EeV assume an inhomogeneous distribution of sources, probably following the large-scale structure distribution in our neighborhood (D. Harari et al. 2014Harari et al. , 2015; P. G. Tinyakov & F. R. Urban 2015; N. Globus & T. Piran 2017; A. di Matteo & P. Tinyakov 2018; C. Ding et al. 2021; D. Allard et al. 2022; T. Bister & G. R. Farrar 2024), or that UHECRs originate from a dominant local source and propagate diffusively through intergalactic magnetic fields (M. Giler et al. 1980;V. S. Berezinskii et al. 1990;S. Mollerach & E. Roulet 2019). The resulting dipole amplitude and direction as a function of energy depend on the source distribution, intervening magnetic fields, and the spectrum and mass composition of the UHECRs. If the sources were distributed like galaxies and if deflections in the Galactic magnetic field could be neglected, a dipolar cosmic-ray anisotropy would be expected in a direction close to that of the dipole associated with the galaxies themselves (P. Erdogdu et al. 2006). On average, a smaller dipolar anisotropy is expected for larger source densities or for more isotropically distributed sources. The amplitude is predicted to grow for increasing energies as a consequence of the shrinking of the horizon due to interactions with the radiation backgrounds (resulting in a larger contribution of the nearby nonhomogeneously distributed sources). Moreover, the Galactic magnetic field modifies the direction and amplitude of the dipolar anisotropy observed at Earth with respect to that of the particles reaching the halo (D. Harari et al. 2010;Pierre Auger Collaboration 2018). These effects depend on the rigidity of the incoming particles and usually lead to a decrease in the dipolar amplitude and also affect the higher multipoles. In particular, for sources following the distribution of galaxies and considering a mean rigidity of the cosmic rays above 8 EeV compatible with the measured mass composition, the direction of the expected dipole would be transformed at Earth after the deflections in the Galactic magnetic field to a direction closer to the observed dipole direction (Pierre Auger Collaboration 2017c).

The energy spectrum and mass composition data measured with the Pierre Auger Observatory above 6 × 10 17 eV can be described by a model with two populations of sources, dominating the flux above and below a few EeV, respectively (Pierre Auger Collaboration 2023a). It has been shown that the data can be well reproduced if the high-energy extragalactic population emits a mixture of heavy and intermediate-mass nuclei with a very hard spectrum and a small rigidity cutoff, while the low-energy component is a mixture of light and intermediate-mass nuclei with a soft spectrum.

In this work, the expected dipole amplitude as a function of energy and the expected directions at Earth are computed for sources emitting cosmic rays according to the best-fitting model for the high-energy population in Pierre Auger Collaboration (2023a) and distributed following the galaxies detected in infrared (IR) in the Two Micron All Sky Survey (2MASS) catalog (M. F. Skrutskie et al. 2006) described in Pierre Auger Collaboration (2022) (with distances from the HyperLEDA database; D. Makarov et al. 2014). The expectation for the dipole amplitude is calculated in the four energy bins above 4 EeV for a population of equal-luminosity sources with number density 10 -5 and 10 -4 Mpc -3 selected randomly from a volume-limited sample of the IR catalog up to 120 Mpc (and from a uniform distribution for larger distances and in the region of the Galactic plane mask). The sources are considered to emit a mixed composition of particles with spectrum E 2 with a rigidity cutoff R cut = 10 18.15 eV (Pierre Auger Collaboration 2023a). The fraction of the total source emissivity above 10 17.8 eV in each mass group is (I H , I He , I N , I Si , I Fe ) = (0., 0.236, 0.721, 0.013, 0.031) (Pierre Auger Collaboration 2023a). The particles Note. For every energy bin, the dipole components d = (d x , d y , d z ), the five independent quadrupolar components Q ij , the quadrupole amplitude Q, and the 99% CL upper limit Q UL are presented. The x-axis is chosen along the α = 0 direction.

were propagated taking into account the interactions with the extragalactic background radiation as described by the Gilmore In Figure 6 the direction of the mean dipole of the simulations and the 68% CL sky region obtained for 10 3 realizations of the source distribution for a density of 10 -4 Mpc -3 are presented.

The results for a source density of 10 -5 Mpc -3 are similar, but the contour regions are larger owing to the larger fluctuations in the source distribution. 101 In the left panel of Figure 7 the median and 68% CL range of the dipole amplitudes are shown for the two source densities considered, also including for For energies below the ankle, the results for the low-energy component of the combined fit of the energy spectrum and mass composition measured with the Pierre Auger Observatory can be described by two different scenarios (Pierre Auger Collaboration 2023a). The first one consists of a Galactic contribution of nitrogen and an extragalactic contribution of protons (which could be produced by interactions of nuclei from the high-energy population in the environment of the sources). The second one consists of an extragalactic contribution of mixed composition (proton, helium, and nitrogen). In both scenarios, there is also a high-energy extragalactic component with mixed composition, as considered before.

The dipolar anisotropy of the nitrogen Galactic component of the first scenario is expected to point close to the Galactic center, in agreement with the measured R.A. phase at low energies, although its amplitude may be close to the present upper bounds at EeV energies (Pierre Auger Collaboration 2012).

Regarding the second scenario, the arrival direction distribution of cosmic rays from the lower-energy extragalactic component is expected to be more isotropic, since at those energies cosmic rays can reach Earth from distances larger than 1 Gpc and the distribution of matter in the Universe at very large scales looks quite uniform. However, even in the case of a perfectly uniform distribution of cosmic rays in some reference frame outside the galaxy, the distribution of arrival directions at Earth is not expected to be isotropic. In particular, the relative motion of the observer with respect to the isotropic reference frame gives rise to an anisotropy that is expected to be dipolar before traversing the Galactic magnetic field, the so-called Compton-Getting effect. If the isotropic reference frame is taken as that of the cosmic microwave background (CMB), this dipolar anisotropy amplitude has been estimated to be of order 6 × 10 -3 (M. Kachelriess & P. D. Serpico 2006), pointing in the direction of the CMB dipole, (l, b) = (264°, 48°) (Planck Collaboration et al. 2020). This amplitude is an order of magnitude smaller than the measured anisotropies at the highest energies and thus represents a subdominant effect above the ankle, but it is of the order of the anisotropies measured below the ankle. The propagation in the Galactic magnetic field distorts the dipolar pattern, changing the amplitude and direction of the dipole observed at Earth, which can have a notable effect on the predictions for the R.A. phase (S. Mollerach et al. 2022). In particular, the dipole direction may shift below a few EeV toward the inner spiral arm direction, at (l, b) ; (-100°, 0°), which can be relevant to explain the observed change in phase. In addition, as a consequence of the Galaxy rotation, the Galactic magnetic field has a small electric component that leads to a directiondependent cosmic-ray acceleration, possibly affecting the anisotropies at a level comparable to the Compton-Getting effect (S. Mollerach et al. 2022).

Clearly, more precise measurements of the anisotropies in the energy range from 0.1 EeV to a few EeV, possibly discriminated by mass composition, and a better knowledge of the Galactic magnetic field will be very useful to shed light on the origin of the cosmic rays below the ankle and on the details of the Galactic-to-extragalactic transition.

Conclusions

The results of the anisotropy searches using the arrival directions of UHECRs detected in Phase 1 of operation of the Pierre Auger Observatory, corresponding to the 19 yr of data gathered before the implementation of the AugerPrime upgrade, were presented. The significance of the established equatorial dipole for the cumulative energy bin above 8 EeV is now 6.8σ, and that for the bin 8-16 EeV is 5.7σ, surpassing the discovery level. For the 8 EeV cumulative energy bin, the statistics are such that we can divide the data set into time-ordered bins, and we find no time variation of the equatorial dipole modulation with an upper limit on the long-term rate of change of 0.003 yr -1 at the 95% CL. If a quadrupolar distribution is allowed in the flux parameterization, the obtained quadrupolar moments are not statistically significant.

The dipole amplitudes increase with energy above 4 EeV. Model predictions were presented for sources following the distribution of galaxies, assuming a source density of 10 -5 Mpc -3 or 10 -4 Mpc -3 , and emitting according to the model for the highenergy population from our combined fit of spectrum and composition (Pierre Auger Collaboration 2023a). The predictions for the dipole amplitude and direction were shown, as well as for the quadrupole amplitude, which are consistent with data within their uncertainties, although some tension with the small observed quadrupole amplitudes seems to be present, in particular for the lowest source density considered and the highest energy bin.

By studying the R.A. distribution of the events, the equatorial component of the dipole down to 0.03 EeV was computed. For energies below 4 EeV the amplitudes are below 1%, compatible with isotropic expectations. However, in most of the energy bins the phase points are consistently close to the Galactic center phase and nearly opposite to the Galactic center phase at energies above 4 EeV. This could be due to the observed anisotropy having a predominant Galactic origin below 1 EeV and a predominant extragalactic origin above a few EeV. Alternatively, this could be caused by the effects of the Galactic magnetic field on an extragalactic flux.

The encouraging prospects of obtaining mass composition estimators on an event-by-event basis with AugerPrime, as well as improved mass estimators with Phase 1 data, will give more information about these anisotropies. An analysis including these estimators to study the different dipole amplitudes 4-8, 8-16, 16-32, 32) EeV.