Covering the celestial sphere at ultra-high energies

Despite deflections by Galactic and extragalactic magnetic fields, the distribution of ultra-high energy cosmic rays (UHECRs) over the celestial sphere remains a most promising observable for the identification of their sources. Thanks to a large number of detected events over the past years, a large-scale anisotropy at energies above 8 EeV has been identified, and there are also indications from the Telescope Array and Pierre Auger Collaborations of deviations from isotropy at intermediate angular scales (about 20 degrees) at the highest energies. In this contribution, we map the flux of UHECRs over the full sky at energies beyond each of two major features in the UHECR spectrum -- the ankle and the flux suppression --, and we derive limits for anisotropy on different angular scales in the two energy regimes. In particular, full-sky coverage enables constraints on low-order multipole moments without assumptions about the strength of higher-order multipoles. Following previous efforts from the two Collaborations, we build full-sky maps accounting for the relative exposure of the arrays and differences in the energy normalizations. The procedure relies on cross-calibrating the UHECR fluxes reconstructed in the declination band around the celestial equator covered by both observatories. We present full-sky maps at energies above 10 EeV and 50 EeV, using the largest datasets shared across UHECR collaborations to date. We report on anisotropy searches exploiting full-sky coverage and discuss possible constraints on the distribution of UHECR sources.


The quest for anisotropies at ultra-high energies
After over a decade of observation with the largest cosmicray observatories at ultra-high energies (UHE, E > 10 18 eV ≡ 1 EeV), the arrival directions of the most energetic particles known to date appear to be nearly isotropic, as measured from the Northern Hemisphere by the Telescope Array (TA) and from the Southern Hemisphere by the Pierre Auger Observatory (Auger). The large exposure cumulated by these two arrays over the years is nonetheless starting to enable a glimpse at preferred directions beyond the ankle of the cosmic-ray spectrum, which may just be the tip of the iceberg. At E Auger > 8 EeV, the Pierre Auger Collaboration reported the detection of a large-angular scale modulation of the cosmic-ray event rate in right ascension, with a small but significant (> 5 σ) amplitude of 4.7 ± 0.8 % * e-mail: biteau(at)ipno.in2p3.fr [1,2]. At energies where the cosmic-ray (CR) flux is suppressed, E TA > 57 EeV, the Telescope Array Collaboration reported an indication (> 3 σ) of clustering of events on an intermediate angular scale of 20 • [3]. The most significant excess found with the Pierre Auger Observatory dataset, for threshold energies between E Auger > 40 EeV and E Auger > 80 EeV, did not reach the 3 σ level using the procedure adopted for model-independent searches [4].
The amplitude and direction of UHE anisotropies is expected to vary with energy. Observations beyond the ankle give access to a non-negligible fraction of the local Universe, e.g. up to 1 Gpc for cosmic rays with E = 10 EeV [5], while the horizon shrinks down to less than 100 Mpc above 100 EeV, being limited by interactions with far-infrared and microwave diffuse photon fields (the CIB and CMB). The highest-energy events possibly arise from a smaller number of sources and should suffer less from deflections within the Galactic magnetic field. 1 The study of anisotropies over the full celestial sphere in these regimes could thus provide invaluable insights into both the propagation of UHECRs over cosmic scales and the distribution of their sources, be they localized in a few directions or distributed following local structures.

Ultra-high energy datasets
Following previous joint searches over the full celestial sphere [7,8], we aim in this contribution at establishing and characterizing datasets covering energy ranges above the ankle and above the flux suppression of the cosmic-ray spectrum. In Sect. 2.1, we describe the datasets collected above threshold energies measured by the Telescope Array and the Pierre Auger Observatory, E TA and E Auger , respectively. We combine in Sect. 2.2 these datasets above common threshold energies, as defined from a match in flux in the declination band covered by both observatories.

Telescope Array and Pierre Auger Observatory datasets
The Telescope Array has been fully operational since May 2008. Data with fiducial cuts described in [9] were shared up to May 2017 above an energy threshold of E TA = 10 EeV, where the array of 507 scintillator detectors is fully efficient. The detectors of the Telescope Array pave an area of nearly 700 km 2 , with a zenith-angle coverage up to 55 • . The geometrical exposure associated with the dataset studied in this contribution, corrected for bin-to-bin migration induced by the limited energy resolution, is estimated to 11,500 km 2 sr yr at E TA > 10 EeV, increasing by less than a percent at the highest energies (11,600 km 2 sr yr for E TA > 50 EeV).
The Pierre Auger Observatory started taking data in 2004 and has been fully operational since January 2008. 1 Median deflections are estimated to be on the order of ∼ 3 • × Z × (E/100 EeV) −1 , with a spread of similar amplitude [6].
Two datasets were shared, with cuts optimized for, on the one hand, analyses above a full-efficiency threshold of E Auger = 4 EeV [1] and, on the other hand, analyses at the highest energies [10].
Following [1], the low-energy dataset consists of events detected from 2004 January 1 up to 2016 August 31. Zenith angles up to 60 • are covered using so-called "vertical" events while a different reconstruction is used from 60 • up to 80 • for so-called "inclined" events. The 1,600 water-Cherenkov detectors of the Pierre Auger Observatory are deployed over area of nearly 3,000 km 2 resulting in a geometrical exposure amounting to 76,800 km 2 sr yr for the period covered by the lowenergy dataset. Correcting for energy resolution effects yields an unfolded exposure of 78,400 km 2 sr yr. At higher energies, above E Auger = 40 EeV [4], the collected dataset covers the period 2004 January 1 up to 2017 April 30 as in [10], with a geometrical exposure of 86,900 km 2 sr yr, that is 91,300 km 2 sr yr after correction for resolution effects. 2

Flux match in the common declination band
The exposures from the Telescope Array and Pierre Auger Collaborations are determined with an accuracy better than 3 %, a level at which the contrast in flux reconstructed from the Northern and Southern hemispheres is not expected to have a significant impact on anisotropy constraints given current UHECR statistics. On the other hand, full-sky anisotropy studies could be substantially impacted by a systematic shift in the energy scales, the uncertainty on which is estimated to be 21 % and 14 % for the Telescope Array and Pierre Auger Observatory, respectively. Shifting the energy scale of a power-law spectrum of index Γ by a factor r results in a shift by a factor r Γ−1 on the flux integrated above a fixed energy threshold. Thus, a 10 − 20 % mismatch in energy scale could result in a 30 − 60 % mismatch in flux for an index Γ ∼ 4. To account for such a flux mismatch, an analysis approach has been devised by the joint Working Group to ensure a proper cross-calibration of the datasets. The approach consists in comparing the flux reconstructed by each instrument in the fraction of the sky covered by both of them. The timeintegrated exposure corresponding to each dataset being practically constant as a function of right ascension over the time period considered, a common declination band can be adopted to perform such a flux estimation. The overlap between the skies covered by the two observatories is illustrated in Fig. 1. Following [8], a common declination band ranging from −12 • to 42 • is used to perform the cross-calibration of the datasets from the Pierre Auger Observatory and Telescope Array.
Because of different variations of the exposure as a function of declination for each observatory, a flux excess localized in the common band could cause different flux estimates. This would particularly be the case for a flux estimator built as d n N( n)/ d n ω( n), where N and ω are the number of events and exposure, respectively, in the direction n. Instead, an unbiased estimator of the flux can be quite naturally obtained as d n N( n)/ω( n) = i 1/ω( n i ), where i indexes the event [8].
The unbiased estimator of the flux reconstructed by each observatory in the common declination band is shown in Fig. 2 as a function of energy threshold. At low energies, driven by the full efficiency of the Telescope Array, we adopt a threshold E TA > 10 EeV. The flux reconstructed from the Telescope-Array dataset above this energy threshold is matched with the Auger dataset for an energy threshold E Auger > 8.86 EeV, as detailed in Table 1. At higher energies, driven by the threshold adopted by the Pierre Auger Collaboration in model-independent anisotropy searches [4] and by indications from model-dependent studies [10], we reconstruct the flux in the common declination band above E Auger > 40 EeV. A compatible flux, reported in Table 1, is obtained from the Telescope-Array dataset above an energy threshold E TA > 53.2 EeV. The amplitude of the mismatch in energy scale shown in Fig. 2 is compatible with that reported in previous efforts, such as [7] where a three times smaller dataset was investigated above E Auger/TA > 8.5/10 EeV, and [8] with a ∼ 40 % smaller dataset above E Auger/TA > 42/57 EeV.
The shift between the energy scales of the two observatories ranges from ± 6 % up to ± 14 %, over the range of interest for this study. The absolute values are within systematic uncertainties reported by each Collaboration. The evolution of the shift with energy, that is the stretch in energy-scale, is estimated to ∼ 10 % per energy decade. Energy-dependent systematic uncertainties have been investigated in the joint Working Group dedicated to the study of the cosmic-ray spectrum, concluding on estimations at the level of ∼ 3 % and < 9 % per energy decade for the Auger and Telescope Array datasets, respectively.

Search for anisotropies at all scales: analysis techniques
Full-sky coverage provides access to an unbiased estimation of the angular power spectrum at all angular scales [7], particularly at low multipoles (dipole, quadrupole) where the full extent of anisotropy patterns can be probed without being impacted by the windowing effect resulting from partial coverage. Beyond a search based on a sphericalharmonics expansion, we also study, for the first time with full-sky coverage, localized excesses at the highest ener-gies. Such searches are made possible by the reconstruction of the UHECR flux and significance maps.

Reconstruction of the flux and angular power spectrum
As discussed in Sect. 2.2 and [8], an unbiased estimator of the flux can be calculated as i 1/ω( n i ). Full-sky maps are built estimating the flux in circular windows centered on a HEALPix 3 grid with nSide = 64. This parameter corresponds to an average pixel size of less than 1 • , which matches the angular resolution of the instruments. Following previous searches [1,8], radii of 45 • and 20 • are adopted as default values for the datasets gathered at E Auger/TA > 8.86/10 EeV and E Auger/TA > 40/53.2 EeV, respectively.
As for the flux, the angular power spectrum can be estimated from the inverse exposure as: where Y m are the spherical harmonics.
Monte-Carlo simulations are used to estimate the distribution of the angular power that is expected to emerge from limited statistics, under the assumption of an isotropic distribution. The spread in the distribution of the angular power can be attributed to two sources: • the limited number of events over the full sphere. This source of uncertainty can be singled out simulating a number of events identical to that gathered over each energy threshold according to the directional exposure presented in Fig. 1; • the uncertainty in the energy cross-calibration procedure, due to the limited number of events in the crosscalibration band. The cumulated effect of this source of uncertainty and of the limited statistics over the full sphere can be estimated by varying, during the Monte Carlo simulation process, the normalization of the exposure of each observatory within an effective uncertainty. The main source of uncertainty is estimated from the resolution on the flux in the common band, as presented in Table 1. 4

Quantifying the deviation from isotropy in all directions
In addition to the flux evaluation, we aim at providing an estimation, over the sphere, of the uncertainty on the flux, whose variance is not only proportional to the flux but also to the exposure, larger in the Southern part of the sky, smaller in the Northern one. More precisely, the significance of deviation from isotropy is estimated in each circular window centered over the HEALPix grid. The Li & Ma estimator [11] (Eq. 17) is used, defining each circular window as an ON-region and the rest of the sky as the OFF-region from which the background rate can be derived. The background rate is normalized by a factor, α, computed as the ratio of the integral exposure over the ON-and OFFregions. By convention, significances corresponding to deficits (N ON < αN OFF ) are quoted as negative values. We checked that the presence of localized anisotropies in the OFF region does not affect in a significant manner the Li & Ma estimator. A comparison with a binomial estimator of the significance, such as used in [4], yields results consistent with the approach followed here. Overall, these tests all converge on an estimation of local significances with an accuracy on the order of 0.2 − 0.3 σ.
The procedure described above only provides estimates of local significances, not accounting for the blind search for flux excesses. The penalization procedure accounting for the search over the HEALPix grid, and over different search radii wherever relevant, is discussed in Sect. 5.

Results with full-sky coverage
The datasets gathered from the Pierre Auger Observatory and the Telescope Array encompass a total of ∼ 31,000 events at E Auger/TA > 8.86/10 EeV and ∼ 1,000 events at E Auger/TA > 40/53.2 EeV, corresponding to an increase by + 200 % and + 60 % with respect to the datasets studied in [7] and [8], respectively.

Mapping the sky beyond the ankle
The flux map reconstructed at E Auger/TA > 8.86/10 EeV is presented in equatorial coordinates in Fig. 3 (left), as estimated in 45 • -radius windows. A large-scale feature is apparent from the enhanced flux in the hemisphere South of the supergalactic plane with respect to the North one. The presence of a smaller scale feature could be inferred from the median flux reconstructed in the North-West quadrant in equatorial coordinates. As presented in Fig. 3 (right), the significance of the North-West feature around (α, δ) ∼ (240 • , 60 • ) is mild with respect to the minimum and maximum significances observed at declinations between ∼ − 60 • and ∼ 0 • . This behavior can be understood from the decrease in exposure with increasing declination.
The angular power spectrum estimated at E Auger/TA > 8.86/10 EeV is shown in Fig. 3, bottom. The largest deviation from isotropy is obtained for the dipole, at = 1, with a local significance of 2.5 σ. As indicated by the dark-and light-gray bands, the uncertainty on the power at = 1 is largely driven by the limited number of events in the common declination band. Such an impact on the dipole can be expected from a varying contrast between the Northern and Southern hemispheres, which impacts the estimation of the component along the Earth rotation axis.

Mapping the sky above the flux suppression
The flux and significance maps reconstructed at E Auger/TA > 40/53.2 EeV using 20 • -radius windows  Fig. 4. Most notably, possible flux enhancements appear along the supergalactic plane, particularly at declinations δ ∼ ± 50 • for the most significant ones.
The angular power spectrum estimated from observations at E Auger/TA > 40/53.2 EeV is shown in Fig. 4, bottom. The largest deviation from isotropy (local 2.8 σ) is obtained for a multipole at = 14, corresponding to an angular scale of about 180 • /l ∼ 13 • . Given that no specific angular scale is preferred a priori at these energies, a penalization for a search over 20 angular scales, or 20 multipoles, results in a post-trial deviation from isotropy at the 1.6 σ level.

Dipolar component beyond the ankle
Constraints on the dipolar component reconstructed at E Auger/TA > 8.86/10 EeV are summarized in Table 2. As also shown in Fig. 3, bottom, no significant deviation from isotropy is obtained from the 3D analysis presented in this work. Nonetheless, the most likely amplitude of the dipolar component along the Earth rotation axis and along the perpendicular plane, as well as associated uncertainties, are provided in Table 2 for the sake of comparison with previous studies. Constraints from Auger data only, presented in the first two lines of Table 2, were derived in [1,2] from a 2D Rayleigh analysis as a function of right ascension and azimuth. Because of partial sky coverage, constraints on the dipole are obtained assuming either no power beyond = 1, or beyond = 2. The larger number of effective free parameters in the 3D analysis followed in this work results in a milder deviation from isotropy than that obtained from a 2D study as a function of right ascension.
A comparison of the central value reconstructed for the component along the Earth rotation axis shows a very good agreement between the 3D and 2D reconstructions. Notably, the limiting factor for the 3D analysis remains the number of events in the cross-calibration band. A somewhat larger tension can be inferred from the estimated central values for the perpendicular component reconstructed in this work and in [1]. This tension is nearly entirely alleviated when allowing for the presence of power in the quadrupolar component at = 2 [2]. Interestingly, the presence of a small-amplitude quadrupolar component could be expected if the sources of UHECRs followed the distribution of local extragalactic matter (e.g. [12]).

Search for overdensities above the flux suppression
Both the Telescope Array and Pierre Auger Collaborations have searched for overdensities in previous studies, with approaches independent of any source model. In These directions are 25 • and 11 • away from the mostsignificant directions obtained in [4] and [3], respectively. Given the difference in energy threshold used in this work and those adopted in [3] and [4], a difference in peak direction on the order of the preferred search radius is considered as a reasonable agreement.
No straightforward comparison with previous searches can be made for the local significance of deviation from isotropy: partial-and full-sky-coverage studies were performed above different energy thresholds and thus exploit quite different number of events given the steep decrease of the UHECR spectrum with increasing energy. As shown in Figs. 5 and 6, the presence of both negative and positive local significances on the order of 4 σ is indicative of a distribution compatible with isotropy. This can also be inferred from a comparison of the distribution of the significance with a Gaussian distribution: both the mean values and spread estimates are compatible at the 1 − 2 σ level with the standard normal distribution expected from isotropy.
We performed a Monte-Carlo computation of the penalty factor to be applied to the local p-values corresponding to the maximum significance reconstructed at 15 • and 20 • : it is estimated to be on the order of 10 4 when accounting for both the scan over the sphere and over the search radius. For comparison, previous search strategies developed independently by each Collaboration would re-    sult in penalty factors on the order of 10 3 and 10 5 for the energy-independent and energy-dependent searches. The two largest reconstructed local significances, penalized for the procedure adopted in this work, correspond to posttrial significances of 2.2 and 1.5 σ.

Conclusion
The Telescope Array and Pierre Auger Collaborations have gathered the largest dataset at ultra-high energies available to date. The combination of data acquired from the Northern and Southern hemispheres enables coverage over the full celestial sphere as demonstrated in previous efforts presented at UHECR conferences. For the first time, studies of both the energy range beyond the ankle and above the flux suppression are presented, where largescale and intermediate-scale anisotropy patterns could be expected.
The number of events gathered at E Auger/TA > 8.86/10 EeV increased by a factor of three with respect to previous effort and by a factor of 1.6 above 40/53.2 EeV. This large dataset enabled a robust cross calibration of the flux in the common declination band, indicating not only a shift between the energy scales of the two observatories but also a stretch, which has been further investigated by the joint Working Group dedicated to the cosmic-ray spectrum. While the full-sky model-independent studies presented in this work did not unveil significant anisotropies, either because of the limited number of events in the common declination band above 8.86/10 EeV or because of the blind-search approach above 40/53.2 EeV, interesting leads could be further pursued. This is particularly the case for constraints on the quadrupolar components in the lowenergy band and for searches along the supergalactic plane at the highest one, where the presence of a "ring of fire" could possibly be inferred from the flux and significance maps shown in Sect. 4.2.
Future searches for anisotropies based on the dataset presented in these Proceedings may shed new light on the distribution of UHECRs over the celestial sphere, be it on large scales beyond the ankle, where the cumulative flux from sources up to 1 Gpc could contribute, or on intermediate angular scales at higher energies, where contributions from a few objects following local extragalactic structures could stand out. The hunt for the most extreme accelerators in the Universe along the anisotropy trail continues.