Measurement of anisotropy and the search for ultra high energy cosmic ray sources

20 eV that are observed through the giant air showers they produce in the atmosphere. These particles carry information on the most extreme phenomena in the Universe. At these energies, even charged particles could be magnetically rigid enough to keep track of, or even point directly to, the original positions of their sources on the sky. The discovery of anisotropy of UHECRs would thus signify the opening of an entirely new window onto the Universe. With the construction and operation of the new generation of cosmic ray experiments—the Pierre Auger Observatory in the Southern hemisphere and the TelescopeArray in the Northern one—the study of these particles, the most energetic ever detected, has experienced a jump in statistics as well as in data quality, allowing for much better sensitivity in searching for their sources. In this review, we summarize the searches for anisotropies and the efforts to identify the sources of UHECRs which have been carried out using these new data.


Introduction
Since the first detection of a cosmic ray (CR) particle with energy in excess of 10 20 eV by J. Linsley at the Volcano Ranch in 1963 [1], the origin of these particles, the most energetic produced in nature, has been an enduring problem in astroparticle physics. These ultra high energy cosmic rays (UHECRs) that come to Earth have a sizeable hadronic cross section: their first interaction in the atmosphere occurs roughly at a column depth of 100 g cm −2 . Given that the proton column density in the Galaxy for most of the UHECR arrival directions is of order 10 −2 g cm −2 , and for the Universe as a whole even smaller (of order 10 −3 g cm −2 Gpc −1 ), the particles arriving at Earth had no prior interactions with protons except perhaps in the source or its vicinity. Adding the requirement of stability and limiting ourselves to known particles, it follows then that the primary particles propagating over cosmological distances should be photons, neutrinos, protons, or nuclei. 1 UHECRs are detected at ground level through the extensive air showers they induce in the atmosphere. Measurements of the PTEP 2017, 12A104 O. Deligny et al. height of the shower maximum exclude to a large extent the presence of photons or neutrinos, so that the bulk of UHECRs detected on Earth are protons and/or nuclei ranging from hydrogen to iron elements. Current experimental progress on the determination of the abundance of each element goes beyond the scope of this review, so that the chemical composition is considered here as possibly ranging from protons to iron nuclei.
There is now evidence, mainly from the stringent upper limits on the UHE photon and neutrino fluxes [3][4][5], that UHECRs are accelerated by electromagnetic processes and are not the products of the decay of supermassive particles as suggested in the top-down scenarios. While the energy spectrum and the chemical composition measurements provide constraints that help in inferring properties of the acceleration processes, identification of the sources can only be achieved by capturing in the arrival directions a pattern suggestive in an evident way of a class of astrophysical objects. This remains, however, a difficult task, mostly because of the very small value of the UHECR intensity (flux per steradian) at Earth-of order 3 × 10 −1 km −2 sr −1 yr −1 above 10 EeV (1 EeV = 10 18 eV)together with the fact that they experience magnetic deflections during propagation. To collect an increased influx of events, a jump in cumulated exposure as well as improved instrumentation have been achieved in the past decade. The Pierre Auger Observatory, located in the province of Mendoza (Argentina) and covering 3000 km 2 , has been allowing since 2004 scrutiny of the UHECR skyexcept in the northernmost quarter-with a total exposure of 6.6 × 10 4 km 2 sr yr [6]. Another scrutiny, mainly of the northern sky, has been provided by the Telescope Array, located in Utah (USA) and covering 700 km 2 , operating since 2008 with a total exposure of 8.7 × 10 3 km 2 sr yr [7]. These latest-generation experiments have allowed an unprecedented sensitivity in searching for anisotropies. The object of this paper is to review the different searches performed in the last decade in the quest to find UHECR sources by making use of the data collected at these observatories.
The general paradigm motivating the search for anisotropies at ultra high energies is presented in detail in Sect. 2. The intervening magnetic fields in extragalactic space and in the Galaxy are responsible for making this task difficult by generically isotropizing the arrival directions. Worse, there are large uncertainties in the estimates of the strength of these fields. The larger deflections are expected to be due to the Galactic magnetic field, whose strength is estimated to be of order several microgauss in the disk of the Galaxy, as inferred from optical and synchrotron polarization measurements and Faraday rotation measures of pulsars and extragalactic sources. Cosmic ray particles with energies in excess of 10 EeV then have a Larmor radius exceeding the dimensions of the Galaxy. Roughly, deflections are expected to behave as 3 • Z(E/100 EeV) −1 , with Z the charge of the CRs. The highest energy particles are thus expected to have sufficient magnetic rigidity to approximately maintain their initial arrival directions provided that the electric charge is not too large. Moreover, the horizon of the highest energy particles ( 60 EeV) is limited as compared to that of particles of lower energies, because the thresholds of interactions with background radiations filling the Universe and leading to large energy losses are then reached. In this way, only the foreground sources are expected to populate the observed sky maps at these energies.
The erasure of the contribution of remote sources provides a natural mechanism to suppress the unresolved isotropic "background," so that UHECRs should originate from the nearby Universe, with source distances not exceeding 100 Mpc or so. With strong nearby sources present, clusters of events could stand out from the isotropic background. Searches for an excess of close pairs of events could reveal the clusters of events on the sky, thus indicating the presence of such sources. Studies of the arrival directions that would be suggestive of such discrete sources without recourse to catalogs of astrophysical objects are summarized in Sect. 3. On the other hand, and even without compelling indications for discrete sources, a correlation between UHECR arrival directions and the positions of a class of astrophysical objects could reveal an anisotropy that would trace the sources. Such a correlation would validate the prospects of studying individual sources with sufficient exposure. Searches for these cross-correlations with different astrophysical objects and, more generally, with the matter distribution in the nearby Universe are presented in Sect. 4.
In addition to the cross-correlation approach, the non-uniform distribution of sources and the corresponding anisotropic distribution of arrival directions of UHECRs may be revealed by studying the decomposition of the observed flux in spherical harmonics. For a sufficient exposure, a dipole moment in the angular distribution could thus be captured, which is due to the density gradient of CRs embedded in the Galaxy. A quadrupole pattern could also arise if nearby sources are distributed along a plane of matter such as the supergalactic plane. Studying the large-scale patterns contained in the arrival direction distribution as a function of energy is thus another important piece of information to decipher. Such studies are reviewed in Sect. 5.
That no prominent source has been detected so far calls into question the understanding of UHECRs prior to the construction of the latest-generation experiments. Magnetic deflections blur the picture to a larger extent than anticipated. This does not preclude that sources may be identified on a collective basis rather than on an individual one in the future, but another jump in statistics appears necessary. Prospects for such charged-particle astronomy are discussed in the conclusion of this review.

UHECR propagation 2.1. Attenuation: The GZK paradigm
Both protons and nuclei of ultra-high energies attenuate when propagating in the intergalactic space due to interactions with the background radiation: cosmic microwave background (CMB), radio, and infrared backgrounds. These interactions occur at the center-of-mass energies accessible to accelerators, so the cross sections are known, or-in the case of nuclear photodisintegration-can be estimated.
Depending on the energy of the primary particle and its nature, different processes dominate the attenuation. For protons, there are two dominant processes: e + e − pair production on the CMB, which is important at energies (1-5) EeV [8], and pion photoproduction on the CMB at energies of around  EeV. At energies around and above 10 EeV, with which we are mostly concerned in this paper, it is thus the second process which is important, causing the so-called Greizen-Zatsepin-Kuzmin (GZK) effect [9,10]. The GZK effect leads to a cutoff in the spectrum at energies higher than 50 EeV; at higher energies it limits the propagation distance for protons to several tens of Mpc.
In the case of nuclei, the dominant energy loss results from nuclear photodisintegration. The nucleon binding energies are of the order of a few MeV with relatively small variation. The photodisintegration process thus becomes efficient when the energy of a background photon boosted in the nucleus rest frame is of the order of several MeV, the energy necessary to split off an individual nucleon, which is the most frequent outcome of the reaction. The relevant parameter governing the reaction is thus the relativistic factor of nuclei.
For light nuclei, the factor is high enough at the energies of interest to induce photodisintegrations on the CMB background, so that the attenuation of light nuclei is much faster than that of protons, qualitatively, the faster the lighter nucleus. For heavy nuclei like iron, the factor remains below 10 10 even at energies as high as 100 EeV. The photodissociation then occurs on infrared background 3  photons with energies about an order of magnitude higher than the typical CMB photon energy. Since these photons are less abundant than the CMB photons, the attenuation of heavy nuclei is relatively slow and is comparable to that of protons. In any case, for both light and heavy nuclei, the result of a photodissociation is most often a nucleon and a lighter nucleus with the same factor, which in turn is subject to further photodisintegration.
The resulting energy loss lengths strongly depend on the energy of the particles. Examples for protons and iron nuclei are shown in Fig. 1, where contributions from the extragalactic background light and the CMB are separated. The cross sections for pair production can be analytically computed via the Bethe-Heitler formula, while those for pion photoproduction have been precisely measured in accelerator-based experiments and can be accurately modeled [21]. In contrast, the cross sections for photodisintegration of nuclei, especially for exclusive channels in which charged fragments are ejected, have only been measured in a few cases so that phenomenological models have to be used to estimate them. A comprehensive estimation of the systematic uncertainties affecting UHECR propagation, mainly due to the photodisintegration cross sections and to the lack of knowledge of the energy spectrum and redshift evolution of the EBL, can be found in Ref. [22].
Overall, the attenuation of UHECRs eliminates the contribution of distant sources to the observed intensity on Earth. Due to this mechanism, all the observed UHECRs should be coming from nearby sources; the higher the energy, the smaller the size of the collection region.

Deflections in cosmic magnetic fields
Trajectories of charged particles are bent by the intervening magnetic fields, so that the flux from a discrete source is spread over a region of the sky, the size of which depends on the rigidity of the particles. The absolute magnitude of deflections is such that a particle of unit charge and energy 100 EeV is deflected by 0.53 • over a distance of 1 kpc in a regular field of magnitude 1 μG, and by 1.8 • over a distance of 50 Mpc in a random field of magnitude 1 nG and correlation length of 1 Mpc.
The extragalactic magnetic fields are assumed to be random and are known quite poorly. From recent measurements of the Faraday rotations of extragalactic sources, the strength of these fields is constrained from above by ∼ 1 nG for coherence lengths larger than 1 Mpc [23,24]. Some simulations indicate that the fields in voids are even smaller. For instance, according to Ref. [25], for protons of 40 EeV propagating over distances of the order of 500 Mpc, the angular deflections should be of 4/26 Downloaded from https://academic.oup.com/ptep/article-abstract/2017/12/12A104/4665684 by guest on 30 July 2018 the order of 1 • except in the directions of galaxy clusters and filaments, where the magnetic fields encountered are more intense and the deflections are thus larger. Hence, given the reduced horizon of particles at the highest energies, the angular deflections in extragalactic space should remain within a few degrees for the majority of nuclei, except for iron having both a larger electric charge and a deeper horizon.
On the other hand, coherent fields with strength up to 3 μG are known to exist in the Galaxy within a disk of 300 pc thickness, roughly described by a structure with arms with reversing field direction between the arms and displaying variations of strength within them [26,27]. Meanwhile, there are uncertainties in the way the field falls off along the direction perpendicular to the disk and in the Galactic halo. The main features are a northerly directed poloidal component falling off slowly with the distance from the disk, and oppositely directed toroidal fields in the halo [26,27]. Depending on the initial direction outside the Galaxy, deflections for protons of 100 EeV are of the order of 2-4 • . Additional turbulent fields with significant variations from arm to arm are also present on correlation lengths of 100 pc. However, since no systematic change in the propagation direction is expected from multiple small deflections induced by these fields, a net root mean square deflection is a few times smaller than the deflections induced by the regular components [28,29].

Searches for clusterings at small and intermediate angular scales
The ultimate goal of CR astronomy is the study of the astrophysical sources producing these particles. As emphasized in the introduction, there is hope for finding discrete sources at the highest energies thanks to the GZK effect, because the isotropic background of CR arrival directions caused by sources distributed throughout the distant Universe is suppressed. With dominating foreground sources in our part of the Universe, clusters of events could be detectable. Searches for clusterings at small and intermediate angular scales performed at the Pierre Auger Observatory and the Telescope Array are presented in this section.

Blind searches for over-densities
To search for over-densities of events over the exposed sky, a simple and popular technique is to build a smoothed sky map by attributing the observed number of events within a circular window with some specific radius to each sampled point on the exposed sky. The probability of the observed number of events in each sample point is then computed from the binomial distribution by estimating the expected number of events for an isotropic distribution within each circular window. A significance sky map is then derived. 2 Conventionally, positive (negative) significances correspond to over-densities (under-densities).
However, by not specifying a priori the targeted regions of the sky where the excesses are searched for as well as the angular window radius and the energy threshold, the probability/significance sky map obtained in this way suffers from the numerous trials performed. In a simple situation in which each trial would be independent from every other, a probability as low as any specified threshold could always be reached by increasing the number of trials. Hence, the number of trials needs to be accounted for in the estimation of the final probability/significance characterizing an excess. In some  sense, this final probability is the original one "penalized" for the various scans performed on the parameters intervening in the analysis. In the toy case aforementioned, the penalization factor would just be the number of trials. In most cases of interest, each trial is not, however, independent of every other. To calculate the penalized probability of an apparent excess of events, Monte Carlo simulations are then the relevant tool allowing for a perfect mimic of the procedure applied to the analyzed data set, reproducing the various correlations between each trial. Mock samples are generated following an isotropic distribution folded into the directional exposure of the experiment considered. The same number of events as in the actual data is generated, keeping the same energy distribution. On each of these mock samples, the set of parameters scanned from the actual data is optimized to capture the most significant excess anywhere on the sampled grid of the exposed sky. The final probability of the original excess (and so its associated significance) is then the number of samples yielding a more significant excess anywhere in the scanned parameter space normalized to the total number of generated samples.
The most up-to-date search for over-densities performed at the Pierre Auger Observatory was reported in Ref. [31], making use of data recorded from January 2004 to March 2014 (corresponding to an exposure of 66, 400 km 2 sr yr). The exposed sky was sampled using circular windows with radii varying from 1 • up to 30 • in 1 • steps, while the energy thresholds were varied from 40 EeV up to 80 EeV in steps of 1 EeV. The resulting (pre-trial) significance sky map is shown in Fig. 2 for energies in excess of 54 EeV in 12 • -radius windows with parameters leading to the maximal significance. The largest departure from isotropy, indicated with a black circle, is characterized by a pre-trial significance of 4.3 σ and is centered at (α, δ) = (198 • , −25 • ), where 14 events are observed against 3.23 expected from isotropy. It is close to the supergalactic plane (shown as the dashed line) and centered at about 18 • from the direction of Centaurus A (shown as the white star). Once penalized for the trials, the probability of this excess is found, however, to be as large as 69%, so that the observed over-density does not provide any statistically significant evidence of anisotropy.
Similar searches have been performed on data collected at the Telescope Array [32]. One notable difference relies on the unique energy threshold used for this analysis, namely 57 EeV, selected from a prior analysis of arrival directions detected at the Pierre Auger Observatory which had initially led to establishing an anisotropy at the 99% confidence level [33], but which was not confirmed with subsequent data [31]. In addition, only oversampling radii of 15, 20, 25, 30, and 35 degrees were used.
Making use of data recorded from May 2008 to May 2013 (referred to as the "five-year data set" hereafter), out of the 72 selected events with zenith angles less than 55 • , an over-density of 19   clustered within a circular window of 20 • radius was observed around the equatorial coordinates (α, δ) = (146.7 • , 43.2 • ) (near the Ursa Major cluster), whereas 4.5 were expected in an isotropic distribution (in other words, 26% of the events were observed to be located within 6% of the sky). The corresponding pre-trial significance for this "hotspot" is 5.1 σ , while the post-trial significance is 3.4 σ .
An updated analysis is now available, using data recorded between May 2008 and May 2015 (the "seven-year data set") [34]. Out of the 109 selected events above 57 EeV and with zenith angles less than 55 • , an over-density of 24 events is observed in a circular window of 20 • around (α, δ) = (148.5 • , 44.6 • ), while 6.88 are expected from isotropic expectations. The position of the excess is centered 1.5 • away from the one found in the previous search. The corresponding pre-trial and post-trial significances for this hotspot are unchanged, namely 5.1 σ and 3.4 σ , respectively. A sky map in equatorial coordinates of the arrival directions of the 109 events with E > 57 EeV is shown in Fig. 3(a). The blue and red points stand for the directions of the UHECRs for the fiveyear and the latest two-year observation periods, respectively. Figure 3 As an alternative approach, the parameters that maximized the original excess can be fixed and used a priori to perform an anisotropy test without penalty factor by making use of the sixth and seventh year data only. In this case, 4 events are observed against 2.31 expected from an isotropic background. The probability of this marginal excess is estimated to be 20%.
Overall, a definite confirmation of the signal has to await additional data with much larger exposure. On the other hand, a survey with full-sky exposure could provide additional information of interest in the quest to find UHECR sources. First attempts to conduct such surveys by means of the metaanalysis of data recorded at the Pierre Auger Observatory and the Telescope Array are currently under way. A mapping of the entire sky for different energy thresholds will be available in the near future.

Autocorrelation function
A clustering of CR events at a certain angular scale might first reveal itself in the autocorrelation function of the events, which measures the cumulative excess of event pairs separated by the given angular scale over the whole field of view, and not necessarily localized around a single reference point as in the approach described in Sect. 3.1. This test is potentially more sensitive than the blind search in a situation where several small excesses of a similar angular size are present: these contribute coherently-get "stacked"-in the autocorrelation function.
The autocorrelation function at a given angular scale ψ can be expressed as (n data − n bg )/n bg , where n data is the number of pairs of data events separated by the angles within the angular bin corresponding to the scale ψ, and n bg is the corresponding number of pairs in the uniformly distributed background, usually calculated by a Monte Carlo simulation. Because of the limited statistics, one usually works not with the correlation function itself, but with the corresponding cumulative quantity F(ψ) = (N data (ψ) − N bg (ψ))/ N bg (ψ), with N data (ψ) and N bg (ψ) being the total number of pairs separated by angles less than ψ in the data and in the simulated background, respectively. If the data and the background are statistically identical, this quantity fluctuates around zero. For large N bg (ψ), positive (negative) values of F(ψ) much larger than one indicate an excess (deficit) of pairs at a corresponding angular scale.
The probability (more precisely, the p-value) of the excess, if that is found, is determined by Monte Carlo simulations by counting how often the number of pairs with angular separation less than ψ in the simulated sets equals or exceeds N data (ψ). In the limit when 1 N bg (ψ) N tot , with N tot being the total number of events in the data, Poisson statistics can be used and F(ψ) gives directly the probability of the excess in equivalent Gaussian sigmas.
In practice, the angular scale ψ is not fixed a priori and is scanned over. Other parameters like the energy threshold in the UHECR data set may also be scanned to maximize the excess. In this case, to derive the global significance of the excess, a penalty factor should be applied to account for the effective number of trials, in a way similar to that described in the previous section.
The autocorrelation function of the events detected at the Pierre Auger Observatory is presented in Fig. 4 (left panel), which shows the p-value (color-coded) as a function of the angular separation ψ and the energy threshold. The energy threshold is scanned in steps of 1 EeV, while the angular separation is scanned in steps of 0.25 • up to 5 • and of 1 • for larger separations. The largest departure from isotropic expectations, leading to a p-value of 0.027, is obtained for a separation angle of 1.5 • and above 42 EeV, where 41 pairs are observed against 30 expected from isotropy (shown by the cross in Fig. 4). Once penalized for the performed scans, it turns out that about 70% of isotropic realizations lead to p-values less than 0.027, so that no self-clustering is captured through this analysis. 8  Similarly, the autocorrelation function of the events detected at the Telescope Array [35] is shown in Fig. 4 (right panel). The p-values are shown as a function of the separation angle ψ for three energy thresholds: E > 10 EeV, E > 40 EeV, and E > 57 EeV. No global significance is given in Ref. [35] for this analysis in view of the small deviation from isotropy. Note, however, that largest deviations from isotropy (the smallest p-values) occur in the set E > 57 EeV at angular scales 20 • -30 • , consistent with the angular size of the hot spot discussed in the previous section.

Searches for correlations with nearby extragalactic matter
The sources of UHECRs being unknown, a plausible hypothesis can still be made about their space distribution: regardless of their nature, they must, at large scales, follow the distribution of baryonic matter. This fact alone may be sufficient to derive the sky distribution of UHECRs as a function of their propagation parameters (composition and magnetic fields, in the first place), which then may be compared to observations to derive constraints on those parameters. The question of sources may thus be disentangled from other unknowns.
At scales of ∼ < 100 Mpc, the matter distribution in the Universe is inhomogeneous. There are large over-densities of matter corresponding to clusters of galaxies, sheets, and filaments, and under-densities corresponding to voids. The three-dimensional positions of the closest of these structures-those within several hundred cubic megaparsecs-are known from complete galaxy catalogs. UHECR sources must trace this distribution to some extent.
A key parameter in this approach is the space density of sources n. In the extreme case that the sources are very rare, a few per (100 Mpc) 3 or less, the source positions will appear random on the sky despite the fact that they follow the matter distribution, because the latter is essentially homogeneous at such large scales. This situation is already disfavored by the existing observations of the very high energy end of the spectrum: in the case where the distance to the closest source exceeds several tens of Mpc, a complete absence of the super-GZK events is expected, or at least a very sharp cutoff.
In the opposite extreme, where the sources are very numerous, they populate the structures proportionally to the total amount of matter. The distribution of sources will thus trace the galaxy distribution, and so should the sky distribution of the UHECR events, after proper accounting for the distance and propagation effects as described in Sect. 2. This is the limit where the UHECR flux, in 9 principle, can be deduced from the galaxy distribution. The nature of the sources would influence the result only slightly through different clustering properties of different types of galaxies-potential acceleration sites of UHECRs.

Cross-correlation analyses
If the deflections of UHECRs are not too big, their arrival directions may show a cross-correlation with the positions of the nearby sources. A search for such a correlation is the most straightforward way to check whether objects of a given class are sources of UHECRs or not.
The correlation function between the UHECR events and a given catalog of objects is calculated by counting the number of pairs (event)-(catalog object) separated by an angular distance within the range defined by a given angular bin, in a way similar to that described in Sect. 3 for the case of the autocorrelation function. As in that case, for reasons of small statistics, one usually considers the cumulative number of pairs, that is, all pairs with separation smaller than that given. The p-value of the excess, if any, is determined either from Monte Carlo simulation in a way similar to that described in Sect. 3, or semi-analytically by calculating numerically the probability that a single UHECR event falls within a given angular distance from any of the sources, and then using the binomial distribution.
As in the case of autocorrelations, if scans over separation angle and/or other parameters are performed, this probability needs to be corrected for the effective number of trials, calculated again by Monte Carlo simulation where synthetic data sets are generated assuming a uniform distribution of the incident particles and passing each set through the same search procedure as the real data. Several analyses of this type have been performed by both Auger and TA collaborations, which we describe below.
An important point to keep in mind is that what is actually tested by any analysis of this type is the hypothesis that the distribution of the UHECR arrival directions is uniform. Low probabilities indicate that this hypothesis is false, but do not imply by itself that the catalog objects are sources of UHECRs: the actual sources may simply trace the distribution of the catalog objects.
The cross-correlation tests with catalogs performed by the Auger collaboration are described in Ref. [31], which we follow here. Several classes of potential sources were considered: the 2MRS catalog of galaxies [36], the Swift BAT catalog of AGNs [37], and a catalog of radio galaxies with jets [38]. For each of the catalogs, the scans were performed over a CR threshold energy of 40 EeV < E th < 80 EeV, an angular scale of 1 • < ψ < 30 • , and a catalog distance cut from 10 to 200 Mpc. In each case the post-trial probability P has been calculated. Figure 5 shows the non-penalized (pre-trial) p-values for the three catalogs considered as indicated on the plots. In each case, the distance cut in the catalog is set to give the minimum p-value; these are cited in the caption. After including the trial factors, the post-trial probabilities equal 8%, 1%, and 1.4% for the 2MRS galaxies, Swift BAT sources, and radio galaxies, respectively. No significant correlations are found.
In the case of the X-ray and radio catalogs, an alternative search was also performed where the maximum distance cut was replaced by the cut on the minimum intrinsic luminosity L. This cut was scanned over. The following post-trial probabilities were found: 1.3% with the cut L > 10 44 erg s −1 for the Swift BAT catalog, and 11% with the cut L > 10 40 erg s −1 for the case of radio galaxies. It was concluded that no significant correlations were observed.
An analogous search has been performed by the TA collaboration with the Northern sky events [39]. The following catalogs of objects were considered: the 3CRR catalog containing radio 10 Fig. 5. Cross-correlation of the UHECR events observed by the Pierre Auger Observatory with the 2MRS catalog of galaxies (left), the Swift BAT catalog of AGNs (middle), and the catalog of radio galaxies with jets (right). Each plot shows the dependence of the p-value (color-coded) on the energy threshold in the UHECR data set and the angular separation. The crosses mark the best (minimum) p-values, which are 1.5 × 10 −3 , 6 × 10 −5 , and 2 × 10 −4 for the three cases from left to right, respectively. galaxies detected at 178 MHz with fluxes greater than 10 Jy [40], excluding the Galactic plane |b| < 10 • ; the 2MRS catalog [36]; the extragalactic subset of the Swift BAT catalog [37] consisting of objects that were detected with a significance greater than 4.8 σ in the energy range of 14-195 keV in the first 58 months of observation by Swift BAT; the compilation of the Swift BAT AGNs detected with at least 5 σ significance in the energy range of 15-55 keV in the first 60 months of observation; the 2LAC set [41] consisting of AGNs detected with at least 4 σ significance in the energy range of 100 MeV-100 GeV in the first 24 months of observation by Fermi-LAT with the exclusion of the Galactic plane |b| < 10 • ; and the VCV catalog [42], which is a compilation of several AGN surveys. In all cases, the maximum redshift cut on the catalog objects, the minimum energy cut on CR events, and the angular scale of correlation were considered free parameters over which the correlation was optimized. The post-trial probability was then calculated by repeating the scanning procedure on the large number of isotropic Monte Carlo sets.
The strongest correlation was found with the Swift BAT AGN catalog, with the energy threshold of 62.2 EeV, angular scale 10 • , and maximum redshift 0.02. The sky map of the TA events and the objects corresponding to these cuts are shown in Fig. 6 The post-trial probability that such a correlation occurs as a fluctuation over the isotropic background was estimated to be 1%, not including the penalty for searching in several catalogs. Thus, no significant correlation with the extragalactic objects was found in the Northern sky either.

Correlations with the large-scale structure
A more elaborate approach is to explicitly take into account the catalog distances and the energy attenuation during propagation. The drawback of this approach, however, is that it involves more unknown parameters. First, the composition of UHECRs is unknown, particularly at the highest energies. While protons and iron nuclei attenuate in a qualitatively similar way, the attenuation of the intermediate nuclei is much faster, as explained in Sect. 2.1. As the attenuation determines the contribution of the remote and thus isotropically distributed sources, the propagation uncertainties affect primarily the overall magnitude of the expected flux variations over the sky.
Second, the magnetic deflections may be large depending on the UHECR composition and energy. These comprise both the regular deflections resulting from the coherent Galactic field, and random ones due to the extragalactic and random Galactic fields (see Sect. 2.2). The latter may be characterized by a single parameter-a typical deflection θ (which, in general, may depend on the direction on the sky); the former require modeling of the coherent Galactic field which, at the moment, involves large uncertainties.
Simplifying assumptions are needed to proceed further. In the TA analysis (see Ref. [43] for details) which we describe now, one assumes a pure proton composition. This is consistent with the TA composition measurements [44], but inconsistent with the Auger composition results [45].
Once a pure proton composition is assumed, the major remaining uncertainty is related to cosmic magnetic fields. A further simplifying assumption is made in the TA analysis that when the magnetic deflections are not too large, they may be accounted for by the random Gaussian smearing of the flux characterized by the single angular scale θ treated as a free parameter. This smearing is supposed to account for the finite angular resolution of the experiment, the random deflections in the Galactic and extragalactic turbulent fields, and the deflections in the regular Galactic field. The latter part of the deflections is not random; however, for small deflections, and for a particular type of statistical analysis that is less sensitive to the coherent nature of deflections than to their overall magnitude, this approach appears a reasonable approximation.
With these assumptions, one can model the expected sky distribution of the UHECR flux by propagating CRs from their sources to the Earth. In the TA analysis [43], the source distribution is assumed to trace the galaxy distribution in the nearby Universe. The latter is obtained from the preliminary version of the 2MRS galaxy catalog [36], which is nearly complete everywhere except in the vicinity of the Galactic plane. The flux-limited subsample with apparent magnitude m ≤ 12.5 is complete out to distances of 250 Mpc; beyond this distance the source distribution is taken as uniform. Each galaxy within 250 Mpc is treated as a UHECR source of fixed intrinsic luminosity and spectrum. Its contribution to the total flux at energies higher than a given threshold E thr is calculated taking into account the distance and attenuation during propagation. Individual contributions are smeared with the two-dimensional Gaussian function of width θ. The weighting correction is made to compensate for non-equal representation of the dim and bright galaxies in a flux-limited sample, as described in detail in Ref. [46]. The fraction of the (uniform) UHECR flux coming from distances larger than 250 Mpc is calculated and added to the total flux.
An example of the flux sky map for the energy threshold E thr = 57 EeV and smearing angle θ = 6 • is shown in Fig. 7   boundaries are chosen in such a way that each band integrates to the same flux. No modulation with the experiment exposure is imposed. One may recognize the known nearby structures; they are described in the figure caption.
Similar maps may be constructed at different energy thresholds and different smearing angles. The effect of changing the smearing angle is obvious. Moving the energy threshold changes the overall contrast of the map, which is encoded in the relative areas occupied by each band. The bands become of equal area in the limit of zero contrast. The higher the energy threshold, the higher the contrast due to suppression of the contributions of remote sources.
In the TA analysis, the statistical significance of the correlation between the expected flux and the actual distribution of events is assessed by means of the "flux sampling" test [43,47]. Given a flux map f (n), with any set of events (real data or Monte Carlo generated) one may associate the set of values of the flux map { f i } = { f (n i )} read off at the positions n i of the events. In this sense the events sample the map, hence the name of the test.
In the case at hand, one wants to know whether the two event sets-e.g. the data and the uniformly generated Monte Carlo set-are distributed in the same way over the sky. If the two event sets are distributed in the same way on the sphere, the two associated sets of flux values also have statistically equivalent distributions: e.g., large (or small) flux values would appear in these two sets equally often. If one finds that this is not the case, the space distributions of the events in the two sets must also be different. The test is binless and does not require that the two sets compared have the same number of events.
In practice, one first generates the expected flux map with given parameters E thr and θ.   For TA this analysis has been performed at three energy thresholds E thr = 10 EeV, 40 EeV, and 57 EeV, and at smearing angles varying from 2 • to 30 • . We show in Fig. 8 the results for E thr = 57 EeV obtained with the 7 years of TA data [35]. For each smearing angle, two hypotheses are tested: that the distribution of events is isotropic, and that the distribution follows the prediction of the LSS model with a given smearing angle. The resulting p-values are represented by the empty squares and filled circles, respectively. The dashed horizontal line marks the 95% C.L. The hypothesis of isotropy is tested many times and is ruled out at about the 3 σ level in most of the individual (statistically dependent) tests. The LSS hypothesis is in fact a different hypothesis at each smearing angle. For most, except maybe the smallest and largest angles, this hypothesis cannot be ruled out by this test (this does not, of course, mean that it is true). The sky map of the TA events in equatorial coordinates, together with the flux expected from the LSS hypothesis superimposed with the TA exposure, are shown in Fig. 9.

Overview of the analysis techniques
Large-angle structures are expected to provide the best anisotropy fingerprints in the arrival direction distributions of UHECR events in the energy ranges in which there are many contributing sources 14 and/or in which the flux from each single source is diffused over a large solid angle due to magnetic deflections. The information on the sources is then contained in the moments of the angular distribution of events, usually expressed in the reciprocal space. Hence, the angular distributions of CRs are generally characterized through the reconstructed moments describing each corresponding angular scale. Prior to reviewing the results, a general reminder on the formalism of the moment reconstruction is provided below.

Harmonic analysis in right ascension
Extensive air shower arrays operate almost uniformly with respect to sidereal time thanks to the rotation of the Earth: the zenith-angle-dependent shower detection is then a function of the declination but not of the right ascension. Thus, the most commonly used technique is analysis in right ascension only, through harmonic analysis of the counting rate within the declination band defined by the detector's field of view [48]. Considered as a function of the right ascension only, the flux of CRs can be decomposed in terms of a harmonic expansion: The customary recipe to extract each harmonic coefficient makes use of the orthogonality of the trigonometric functions. Modeling any observed arrival direction distribution, (α), as a sum of N Dirac functions over the circle, (α) = i δ(α, α i ), the coefficients can be estimated from the discrete sums: Here, the recalibrated harmonic coefficients a c n ≡ a c n /a 0 and a s n ≡ a s n /a 0 are directly considered, as is traditionally the case in measuring relative anisotropies. Over-lined symbols are used to indicate the estimator of any quantity. The statistical properties of the estimators {a c n , a s n } can be derived from the Poissonian nature of the sampling of N points over the circle distributed according to the underlying angular distribution (α). In the case of small anisotropies (i.e., |a c n /a 0 | 1 and |a s n /a 0 | 1), the harmonic coefficients are recovered with an uncertainty such that σ c n (a c n ) = σ s n (a s n ) = √ 2/N . For an isotropic realization, a c n and a s n are random variables whose joint probability density function (p.d.f.), p A c n ,A s n , can be factorized in the limit of a large number of events in terms of two Gaussian distributions whose variances are thus σ 2 = 2/N . For any n, the joint p.d.f. of the estimated amplitude, r n = (a c2 n +a s2 n ) 1/2 , and phase, φ n = arctan (a s n /a c n ), is then obtained through the Jacobian transformation: From this expression, it is straightforward to recover the Rayleigh distribution for the p.d.f. of the amplitude, p R n , and the uniform distribution between 0 and 2π for the p.d.f. of the phase, p n . Overall, this formalism provides the amplitude of each harmonic, the corresponding phase (right ascension of the maximum intensity), and the probability of detecting a signal due to fluctuations of an isotropic distribution with an amplitude equal to or larger than the observed one as P(> r n ) = exp (−N r 2 n /4). The first harmonic amplitude and phase, corresponding to the case n = 1, are of special interest 15 and generally draw the particular attention of observers, since they constitute generic expectations from various models. We will discuss the results for this harmonic obtained for EeV and trans-EeV energies when presenting Fig. 10.
Note that the aforementioned formalism can be applied off the shelf only in the case of an exposure that is purely uniform in the right ascension, a condition that is generally not fulfilled. At the sidereal time scale, the directional exposure of most observatories operating with high duty cycle (e.g., surface detector arrays) is, however, only moderately non-uniform. Different approaches are then adopted in the literature to account for the non-uniformities. Defining ω(α) as the directional exposure integrated in declination, a widely used and simple recipe is to transform the observed angular distribution (α) into the one that would have been observed with a uniform directional exposure, (α)/ω r (α), with ω r the dimensionless relative directional exposure defined as ω r (α) = 2πω(α)/ , with the total exposure. In that way, the discrete summations in Eq. 2 are changed as follows: The uncertainty on the recovered coefficients then reads, still in the case of small anisotropies [49], as For variations of ω r (α) of a few percent, these expressions are accurately approximated by σ c n = σ s n = 2/Ñ , so that in such cases, the p.d.f. of the amplitude remains a Rayleigh distribution with the parameter 2/Ñ , while the probability of detecting a signal due to fluctuations of an isotropic distribution with an amplitude equal to or larger than the observed one is P(> r n ) = exp (−Ñ r 2 n /4).

Multipole expansion in right ascension and declination
In general, and in contrast to the simplified approach presented in the last subsection, the flux of CRs (n) can depend on both the right ascension and the declination and thus be decomposed in terms of a multipolar expansion in the spherical harmonics Y m (n): Non-zero amplitudes in the modes arise from variations of the flux on an angular scale 1/ radians. With full-sky but non-uniform coverage, the customary recipe for decoupling directional exposure effects from anisotropy ones consists in defining the recovered coefficients as [50] a m = 4π dn ω(n) with ω(n) the directional exposure providing the time-integrated surface of the experiment to each direction of the sky, and dN (n)/dn the observed angular distribution. Modeling this latter distribution as a sum of Dirac functions in the directions of each event, dN (n)/dn = i δ(n, n i ), the coefficients 16 can be estimated in practice asā where ω r (n) is the relative directional exposure function normalized here to unity at its maximum, and f 1 = dn ω r /4π is the covered fraction of the sky. In the case of small anisotropies, the uncertainty σ m on each a m multipole reflects the Poisson fluctuations induced by the finite number of events: However, with ground-based observatories, coverage of the full sky is not presently possible with a single experiment. The partial-sky coverage of ground-based observatories prevents the multipolar moments a m being recovered in the direct way just presented. This is because the solid angle on the sky where the exposure is zero prevents one from making use of the completeness relation of the spherical harmonics. Indirect procedures have to be used, one of them consisting in considering first the "pseudo"-multipolar moments and then the system of linear equations relating these pseudo-moments to the real ones: Formally, the coefficients a m appear related toã m through a convolution. 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 function ω(n). Assuming a bound max beyond which a m = 0, these relations can be inverted allowing recovery of the moments a m . However, the obtained uncertainty on each moment does not behave as expressed in Eq. 9. In contrast, the uncertainty σ m on the recovered a m coefficients behaves as [51] The inverse matrix K −1 is indexed here by the bound max , because of the dependence of the matrix coefficients on this parameter. In numbers, it turns out that σ m increases exponentially with max . This dependence is nothing but the mathematical translation of it being impossible to know the angular distribution of CRs in the uncovered region of the sky. In most of the practical cases reviewed in the next subsections, the small values of the energy-dependent a m coefficients combined with the available statistics in the different energy ranges do not allow for an estimation of the individual coefficients with a relevant resolution as soon as max > 2. From the estimation of the spherical harmonic coefficients, a more geometric and more intuitive representation of the dipole and quadrupole moments is generally used: 1 + r d · n + λ + (q + · n) 2 + λ 0 (q 0 · n) 2 + λ − (q − · n) 2 + · · · .  [52]. Amplitudes are also reported (black squares) in the two energy bins when the corresponding p−value expected from isotropy is below 10 −3 . Right: Corresponding phases as obtained at the Pierre Auger Observatory (filled black squares) [52], at the Telescope Array (open red circles) [53], and at the Yakutsk array (open blue squares) [54,55].
In this picture, the dipole moment is thus characterized by a vector, whose amplitude r and two angles of the unit vector d are related to the a 1m coefficients through 3 r = √ 3 a 00 a 2 10 + a 2 11 + a 2 1−1 The amplitude r corresponds to the anisotropy contrast of a dipolar flux. The quadrupole, on the other hand, is characterized by the amplitudes (λ + , λ 0 , λ − ) and unit vectors (q + , q 0 , q − ) which are the eigenvalues and eigenvectors of a second order, traceless, and symmetric tensor Q whose five independent components are related to the a 2m coefficients through The eigenvalues are ranked from largest to smallest and assigned to the vectors (q + , q 0 , q − ) that form the principal axes of the coordinate system. The traceless condition of the quadrupole tensor Q forces the relation λ + + λ 0 + λ − = 0 to be satisfied, so that only two of these amplitudes are independent. Hence, two diagnostic parameters are used to characterize a quadrupole anisotropy: the quadrupole magnitude that takes on the value λ + , and the anisotropy contrast of a quadrupolar flux β = (λ + − λ − )/(2 + λ + + λ − ). The orientation of the quadrupole is then described by the three Euler angles determined from the eigenvectors corresponding to each of the principal axes and characterizing the orientation of these principal axes with respect to some reference coordinate system.

Harmonic/multipolar analyses from single experiments
Scrutiny of the large-scale distribution of arrival directions of UHECRs provides important information in the EeV energy range. A time-honored picture is that the ankle is a feature in the energy spectrum that marks the transition between Galactic and extragalactic CRs [56]. In the energy range just below the ankle energy, the propagation regime between the diffusive particle transport and the deterministic flow of particles is expected to induce large-scale anisotropies shaped by the distribution of sources in the disk of the Galaxy and the structure of the coherent Galactic magnetic field. On the other hand, the eventual anisotropies of EeV-scale extragalactic CRs should not reflect the geometry of the Galaxy. These generic benchmark scenarios provide some signatures that should help in establishing the energy at which the flux of extragalactic CRs starts to dominate the energy spectrum. Answering this long-standing question would constitute an important step forward toward understanding the origin of UHECRs. Measurements of the amplitudes r 1 and phases φ 1 of the first harmonic in the right ascension at EeV and trans-EeV energies are shown in Fig. 10. The strongest constraints on the amplitudes are currently provided by the Pierre Auger Observatory [52]. As in none of the energy bins are the p-values for the amplitudes at the level of discovery, the upper limits, at the 99% C.L., are shown in the left panel. Amplitudes are also shown in the two energy bins where the p-values are 1.5 × 10 −4 (1 < E/EeV ≤ 2) and 6.4 × 10 −5 (E > 8 EeV) [58]. Together with these relatively low p-values, an apparent consistency in the phases as measured at the Pierre Auger Observatory [52], at the Telescope Array [53], and at the Yakutsk experiment from an older analysis [54,55] is observed, even though the significances of the corresponding amplitudes are relatively small. Note that the uncertainties on the phases are estimated as √ 2/N /r 1 , so that these uncertainties do not reflect the cumulated statistics only. As already pointed out by Linsley a long time ago, the observed consistency is potentially indicative of a real underlying anisotropy, because consistency of the phase measurements in ordered energy intervals is indeed expected to be revealed with a smaller number of events than needed to detect the amplitude with high statistical significance [59]. Interestingly, the phases derived from experiments operating mainly in the PeV-EeV energy range show a consistent tendency to align in the right ascension of the Galactic center. In this regard, the change of phase observed above 1 EeV toward, roughly, the opposite of the one at energies below 1 EeV provides interesting information.
To further characterize the EeV and trans-EeV angular distributions of CRs, a thorough search for large-scale anisotropies in terms of dipole and quadrupole moments was conducted by the Pierre Auger Collaboration [60]. Assuming that the anisotropic component of the angular distributions is dominated by pure dipole or dipole and quadrupole patterns, searches for significant moments were performed. Within the statistical uncertainties, no strong evidence of any significant amplitude could be captured. For a pure dipolar flux, the reconstructed directions are shown in orthographic projection in Fig. 11 with the associated uncertainties, as a function of the energy. The same change of phase in right ascension as in the case of the first harmonic analysis is observed. In addition, the reconstructed declinations are observed to be in the equatorial southern hemisphere.
The obtained upper limits on dipole and quadrupole amplitudes, derived at 99% C.L., are shown in Fig. 12. The dipole amplitude plotted in the left panel corresponds to the quantity r defined in Eq. 14, characterizing the anisotropy contrast of dipolar-shaped flux. The quadrupole amplitude λ + plotted in the right panel provides the magnitude of similar excesses at antipodal points in the sky. It corresponds to the greatest eigenvalue of the quadrupole tensor presented in Eq. 15. The bounds on the dipole amplitudes as a function of energy are shown in the left panel along with generic estimates of the dipole amplitudes expected from stationary Galactic sources distributed in the disk considering two extreme cases of single primaries: protons and iron nuclei. Both the strength and the structure of the magnetic field in the Galaxy, known only approximately, play a crucial role in the propagation of  Fig. 11. Reconstructed declination and right ascension of the dipole moment at the Pierre Auger Observatory, as a function of the energy, in orthographic projection [57].  [57]. Some generic anisotropy expectations from stationary Galactic sources distributed in the disk are also shown, for two distinct assumptions on the CR composition.
CRs and thus in the predictions of their arrival directions. The best up-to-date models can be found in Refs. [26,27]. As an illustrative case, the bisymmetric spiral structure model with anti-symmetric halo with respect to the Galactic plane described in Ref. [26] was considered in this study, on top of a turbulent field generated according to a Kolmogorov power spectrum. This example is an illustration of the potential power of these observational limits on the dipole anisotropy to exclude the hypothesis that a light component of EeV-CRs comes from stationary sources densely distributed in the Galactic disk and emitting in all directions. Furthermore, assuming that the angular distribution is modulated by a dipole and a quadrupole, the 99% C.L. upper bounds on the quadrupole amplitude corresponding to the principal axis excess that could result from fluctuations of an isotropic distribution are shown in the right panel together with expectations considering the same astrophysical scenario. The same conclusion holds. Hence, the percent limits on the amplitude of the anisotropy exclude the presence of a large fraction of Galactic protons at EeV energies. Similar conclusions were recently obtained with data recorded at the Telescope Array by fitting the angular distribution observed at EeV energies to benchmark sky maps obtained by adding the flux expected from stationary Galactic sources of protons on top of a fraction of isotropic flux [61]. The exact modeling for the expected flux of Galactic protons, taken from Ref. [62], is similar to that described above. An example of a benchmark sky map for energies between 1 and 3 EeV is shown in the left panel of Fig. 13, while the sky map observed in this energy range is shown in the right panel. At most, not more than 1% of the observed flux can be described by the modeled sky maps [61].
Accounting for the inference from the depth of air shower maximum data from both the Pierre Auger Observatory and the Telescope Array that protons are in fact abundant at those energies [44,63], the lack of strong anisotropies provides some indication that this component of protons is extragalactic, gradually taking over a Galactic one. The low level of anisotropy could then be the result of the addition of two vectors with opposite directions, naturally reducing the amplitudes and producing the change of phase observed at EeV energies. This scenario is to be explored with additional data. Increased statistics is thus necessary to probe the anisotropy contrast levels that may exist in this energy range and contain valuable information about the long-standing question on the way the transition between Galactic and extragalactic CRs occurs. Also, a current limitation of the measurements is that neither spectra nor anisotropies can yet be studied as a function of the mass of the particles with adequate statistical precision, measurements that would allow a distinction between Galactic and extragalactic angular distributions.

Multipolar analysis above 10 EeV with full-sky coverage
Above 10 EeV, the whole flux of UHECRs is expected to be of extragalactic origin. Although the actual sources of UHECRs are still to be identified, their distribution in the sky is expected to follow, to some extent and as already stressed in previous sections, the large-scale structure of the matter in the Universe. It is thus interesting to highlight again that above 8 EeV, the amplitude shown in Fig. 10 has a p-value as low as 6.4 × 10 −5 . Assuming that the only significant contribution to the anisotropy is from a dipolar pattern, the amplitude of this signal converts into a (7.3 ± 1.5)% dipole amplitude [58]. This hint may constitute in the near future the first detectable signature of extragalactic CRs observed on Earth.
To characterize further the angular distribution above 10 EeV, the dipole moment on the sphere is of special interest. An unambiguous measurement of this moment as well as of the full set of spherical harmonic coefficients requires full-sky coverage. Currently, this can be achieved piecemeal by combining data from observatories located in both the northern and southern hemispheres. To this 21  end, a joint analysis using data recorded at the Pierre Auger Observatory and the Telescope Array above 10 EeV has been performed in Refs. [64,65]. Thanks to the full-sky coverage, the measurement of the dipole moment reported in these studies does not rely on any assumption on the underlying flux of CRs.
The main challenge in combining the data sets is to account adequately for the relative exposures of both experiments. A band of declinations around the equatorial plane is exposed to the fields of view of both experiments, namely for declinations between −15 • and 25 • . This overlapping region has been used to design an empirical procedure to get a relevant estimate of the relative exposures: for an isotropic flux, the integrated energy spectra measured independently by both experiments in the common band would have to be identical. The commonly covered declination band could thus be used for cross-calibrating empirically the energy spectra of the experiments and for delivering an overall estimate of the relative exposures. Since the shapes of the exposure functions are not identical in the overlapping region, the observed energy spectra are not expected to be identical in the case of anisotropies. For small anisotropies, however, this guiding idea can nevertheless be implemented in an iterative algorithm finally delivering estimates of the relative exposures and of the multipole coefficients at the same time. The uncertainties on the recovered coefficients, however, are larger than expected from Eq. (12) due to the effect of the uncertainty in the relative exposures of the two experiments. This propagation of uncertainty mainly impacts the resolution in the dipole coefficient related to variations of the flux in declination.
The resulting entire mapping of the celestial sphere has revealed a dipole moment with an amplitude r = (6.5 ± 1.9)%, captured with a chance probability of 5 × 10 −3 [65]. No other deviation from isotropy has been observed at smaller angular scales. The recovered moment is visualized in Fig. 14, where the average flux smoothed out at an angular scale of 60 • per solid angle unit is displayed using the Mollweide projection, in units of km −2 yr −1 sr −1 . This map is drawn in equatorial coordinates. The direction of the reconstructed dipole is shown as the white star.
Large-scale anisotropies of CRs with energies in excess of 10 EeV are closely connected to the sources and the propagation mode of extragalactic UHECRs, see e.g., Refs. [66,67]. Due to scattering 22 in the extragalactic magnetic fields, large deflections are expected even at such high energies for field amplitudes in the nanogauss range and extended over coherence lengths of the order of one megaparsec, or even for lower amplitudes if the electric charge of UHECRs is large. For sources distributed in a way similar to the matter in the Universe, the angular distribution of UHECRs is then expected to be influenced by the contribution of nearby sources, so that the Milky Way should be embedded into a density gradient of CRs that should lead to at least a dipole moment. The contribution of nearby sources is even expected to become dominant as the energy of CRs increases due to the reduction of the horizon of UHECRs induced by energy losses that are more important at higher energies. Once folded through the Galactic magnetic field, the dipole pattern expected from this mechanism is transformed into a more complex structure presumably described by a lower dipole amplitude and higher-order multipoles. However, in these scenarios, the dipole moment could remain the only one within reach given the sensitivity of the current generation of experiments. On the other hand, the detection of significant multipole moments beyond the dipole could be suggestive of non-diffusive propagation of UHECRs from sources distributed in a non-isotropic way.

Conclusions and outlook
At ultra-high energies, the popular idea that cosmic-ray astronomy may be feasible relies on two assumptions: on the possibility that the messengers have high enough rigidities to beat the magnetic blurring, and that the reduced particle horizons caused by energy losses eliminate or suppress the isotropic background from distant sources. However, for the same reason, the pattern recognition of the astrophysical sites harboring the accelerators which could be performed from UHECR arrival directions is made difficult by the very small intensity of these particles. It is also possible that the rigidity of the particles is not sufficient: the properties of intervening magnetic fields are uncertain, and the mass composition of UHECRs in the highest-energy region remains largely unknown.
Thorough searches for various possible anisotropies of UHECRs conducted by the Pierre Auger and Telescope Array collaborations during the past decade have been summarized in this review. The data show a remarkable degree of isotropy at all energies, with only a few hints on a possible anisotropic distribution, like large-scale patterns in the EeV energy range, and the concentration of events (the hot spot) in the Northern hemisphere at the highest energies E > 57 EeV (TA energy scale). The large-scale pattern above 8 EeV (Auger energy scale) is best revealed through the first harmonic analysis in right ascension of the Auger data. These indications still require confirmation with larger statistics.
With the data above 55 EeV released by both collaborations, it is worth mentioning that other explorations of possible indications of anisotropies were performed by numerous authors outside of these collaborations which were not covered in this review, e.g. searches for clustering of events [68][69][70]; searches for correlations with extragalactic matter [71][72][73][74][75][76][77][77][78][79]; studies of the lensing effect of the Galactic magnetic field to probe whether the observed arrival directions could correlate with extragalactic objects once unfolded [80,81]; and others. Overall, although a few intriguing correlations with positions in the sky tracing special astrophysical environments have been uncovered, none of the corresponding significances is large enough to provide a compelling signal of anisotropy at present, especially in view of the multiple searches conducted.
The apparent isotropy can be explained by several factors. One of them is the possibility that UHECRs are composed of a mix of light and intermediate/heavy nuclei. In addition, the magnetic field in the halo of the Galaxy may be underestimated due to the low density of electrons in these regions, so that even for protons, the magnetic deflections could still be large at the highest energies. 23 The extragalactic fields, though bound by ∼nG at large scales, may be larger in the ∼Mpc vicinity of our Galaxy in the case that it is embedded in a filament.
The absence of obvious bright sources on the UHECR sky, if not due to large magnetic deflections, may indicate that the sources are too numerous and individually weak. This can be quantified by setting lower bounds on the density of sources. Roughly, the absence of repeaters implies that the number of contributing sources has to be larger than the square of the number of events [82]. Using the events above 70 EeV detected at the Auger Observatory, the best up-to-date derived bounds, whose validity domains are limited for magnetic deflections smaller than the angular scale ψ used to search for the repeaters, range from 7 × 10 −4 Mpc −3 at ψ = 3 • up to 2 × 10 −5 Mpc −3 at ψ = 30 • [83]. For comparison, the density of galaxies is roughly 10 −2 Mpc −3 .
First attempts to get full-sky surveys at ultra-high energies are being developed between the two collaborations. Full-sky coverage is obviously advantageous to probe all possible sources of UHECRs, and is indispensable for the harmonic analysis aimed at revealing possible anisotropies at large angular scales. Also, a multi-messenger approach could give clues for deciphering the origin of the cosmic ray particles. In this context, a full-sky study conducted in a collaboration between Auger, Telescope Array, and IceCube has recently reported on the search for correlations between UHECRs and very high energy neutrino candidates detected by IceCube [84]. It is interesting that the smallest post-trial p-values (corresponding to significances slightly greater than 3 σ ) are obtained when considering the correlations between the directions of cascade events observed by IceCube and those of the UHECRs on an angular scale of 20 • . With increased statistics, this kind of meta-analysis will help us to understand whether or not a contribution to the neutrino signal observed by IceCube arises from the sources of the observed UHECRs.
Present detector exposures are already providing us with important constraints on the origin of UHECRs. Even larger exposures will lead to much better constraints and hopefully will eventually make possible the detection of the brightest sources. The Telescope Array collaboration is building an extension of the surface detector array, reaching a detection surface of about 3000 km 2 . This will boost the UHECR statistics in the Northern hemisphere which, apart from the sensitivity to sources in the Northern sky, is crucial for the all-sky surveys. At the Auger Observatory, upgraded instrumentation is being deployed to equip the detectors with an additional plane of 4 m 2 of plastic scintillators above each station. This will provide us with additional high-statistics measurements of the showers, helping to produce composition-sensitive observables in an energy range including the highest energies, and will open the possibility of performing mass-discriminated anisotropy searches if the composition is mixed with light and heavy elements [85]. The operation of such an upgraded array is anticipated between 2018 and 2024.