Issue 
A&A
Volume 652, August 2021



Article Number  A41  
Number of page(s)  11  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202038459  
Published online  06 August 2021 
Detecting ultrahighenergy cosmic ray anisotropies through harmonic crosscorrelations
^{1}
CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182 21 Prague, Czech Republic
email: federico.urban@fzu.cz
^{2}
Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy
^{3}
INFN – Istituto Nazionale di Fisica Nucleare, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
^{4}
Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
Received:
21
May
2020
Accepted:
2
May
2021
We propose an observable for ultrahighenergy cosmic ray (UHECR) physics: the harmonicspace crosscorrelation power spectrum between the arrival directions of UHECRs and the largescale cosmic structure mapped by galaxies. This crosscorrelation has not yet been considered in the literature, and it permits a direct theoretical modelling of the main astrophysical components. We describe the expected form of the crosscorrelation and show how, if the distribution of UHECR sources traces the largescale cosmic structure, it could be easier to detect with current data than the UHECR autocorrelation. Moreover, the crosscorrelation is more sensitive to UHECR anisotropies on smaller angular scales, it is more robust to systematic uncertainties, and it could be used to determine the redshift distribution of UHECR sources, making it a valuable tool for determining their origins and properties.
Key words: astroparticle physics / largescale structure of Universe
© ESO 2021
1. Introduction
Ultrahighenergy cosmic rays (UHECRs), which impact the atmosphere of the Earth with energies in excess of 1 EeV (10^{18} eV), have remained a mystery since their discovery 59 yr ago (Linsley et al. 1961; Alves Batista et al. 2019). We do not know what they are: Observational data cannot yet fully distinguish between several variants of pure and mixed primary compositions (Castellina 2020; Bergman 2019). We do not know where they come from: The astrophysical sources that generate and accelerate UHECRs have not yet been identified, nor has the type of acceleration mechanism responsible for their formidable energies been discovered (Kotera & Olinto 2011).
What we do know is that the highestenergy rays are most likely extraGalactic. First, if UHECRs were produced within the Galaxy, their arrival directions in the sky would be very different from what we observe (Tinyakov et al. 2016; Abbasi et al. 2017; Aab et al. 2017). Second, barring a cosmic conspiracy that puts an end to the injection spectrum at that very energy, UHECR interactions with cosmological background photons produce a sharp cutoff (the GreisenZatsepinKuzmin limit) in the spectrum corresponding to ∼60 EeV (Greisen 1966; Zatsepin & Kuzmin 1966), and a cutoff is indeed observed in the data (Abbasi et al. 2008; Abraham et al. 2008).
If the sources of UHECRs are extraGalactic, they most probably correlate with the largescale distribution of matter (largescale structure, or LSS). The interactions with the background cold photons limit UHECR propagation to roughly a few hundred megaparsecs (for a review, see Kotera & Olinto 2011). Therefore, the UHECR flux distribution in the sky should be to some extent anisotropic since below 100 Mpc, roughly comparable with the scale of homogeneity expected in the standard cosmological model, the LSS is anisotropic (Pan & Coles 2000; Scrimgeour et al. 2012; Alonso et al. 2015).
How the anisotropy of UHECR sources manifests itself in the observed flux on Earth therefore depends on the original anisotropy of the sources, the UHECR chemical composition, and the properties of intervening magnetic fields – Galactic (GMF) and extraGalactic (xGMF) – that deflect UHECRs and distort the original anisotropic patterns. Chemical composition and magnetic fields are degenerate when it comes to UHECR deflections since such deflections depend on ZB/E, where Z is the atomic number, B is the strength of the magnetic field, and E is the UHECR energy: doubling the field strength is equivalent to doubling the charge (or halving the energy). Chemical composition instead is the only factor that determines the UHECR propagation length at a given energy: different nuclei come from different portions of the Universe and carry different anisotropic imprints, but the relationship between the two is nonmonotonic and nontrivial (see e.g., d’Orfeuil et al. 2014; di Matteo & Tinyakov 2018).
To a large extent, the statistics of the anisotropies in the distribution of UHECRs can be characterised by the UHECR angular autocorrelation (AC), which, in harmonic space, takes the form of the angular power spectrum coefficients C_{ℓ}. Here, the ℓth multipole quantifies the variance of the anisotropies on angular scales θ ∼ π/ℓ (Sommers 2001; Tinyakov & Urban 2015, see our Appendix A for further details). To date, the number of UHECRs collected at the highest energies is low – of the order of a hundred above the cutoff (Alves Batista et al. 2019). Because of this, the UHECR flux is dominated by Poisson statistics: the AC is mostly determined by shot noise, making the underlying correlation with the LSS very hard to detect. Indeed, the indications for anisotropy in the data are tenuous: save for a lowenergy dipole (Aab et al. 2017) and a highenergy hot spot (Abbasi et al. 2014), the angular distribution of UHECR arrival directions appears to be nearly isotropic (di Matteo et al. 2020). Moreover, no anisotropies have been detected at small scales, ℓ ≳ 10 (Deligny 2015; di Matteo et al. 2018), although there are hints at intermediate scales (Caccianiga 2020).
In this work we quantify the possibility of detecting the anisotropy in the UHECR flux through the harmonicspace power spectrum of the crosscorrelation (XC) between UHECR counts and the distribution of galaxies. Such a XC technique was previously proposed to study the anisotropy of the γray sky by Camera et al. (2013), Fornengo & Regis (2014), and Pinetti et al. (2019) and proved successful for several tracers of the LSS (Fornengo et al. 2015; Cuoco et al. 2015; Branchini et al. 2017; Ammazzalorso et al. 2020). A search for a XC between UHECRs and highenergy photons was performed in Alvarez et al. (2016). If UHECR sources statistically trace the LSS, then the positions of these sources, and the arrival directions of UHECRs (if not strongly affected by intervening magnetic fields), should have a nonzero correlation with a galaxy sample up to a given distance. Therefore, the detection or nondetection of the XC signal with galaxies at different redshifts would allow us to test whether UHECR sources are distributed according to the LSS, and to quantify the extent to which the UHECR transfer function, determined by energy losses and intervening magnetic fields, is independent from direction.
There are at least three features that differentiate the XC from other methods (see for instance Koers & Tinyakov 2009 and references therein). First, systematic uncertainties of different ‘messengers’, or observables, should not crosscorrelate, and, under some conditions, statistical noise should not strongly crosscorrelate either. This is because different experiments are different machines exploiting different physical effects. However, within a single dataset, for instance the set of arrival directions of UHECRs, the AC of the noise and systematic errors for that set are certainly nonzero and contribute to hiding any underlying ‘true’ signal. Examples of these systematics for UHECRs would be perturbations in the arrival directions due to deflection by the GMF^{1} or spatial fluctuations in the energy calibration that give rise to leakage from different energy bins^{2}. Thus, in this sense the XC is an experimentally cleaner observable.
Second, in the limit where the UHECR sources are numerous but UHECR detections themselves are not, we can assume that we observe at most one UHECR per source (as seems to be the case given the lack of obvious UHECR multiplets; Abreu et al. 2012; Aab et al. 2019). The much higher number of galaxies leads to a significant improvement in the signaltonoise ratio (S/N) of this XC (see the discussion in Sect. 3). This effectively allows us to probe the anisotropies on smaller scales through the XC than through the AC, underlying the importance of using both observables.
There are several reasons why those smaller scales (ℓ > 10) are interesting. The experimental angular resolution of UHECR events is around 1°, which corresponds to ℓ ∼ 200: From an experimental perspective, we are not fully taking advantage of the data we already have. Furthermore, smallscale power in the LSS angular distribution is comparable to that at large scales: if UHECRs bear the imprint of the LSS, this smallscale power is not completely suppressed by the GMF, especially once the structured component of the GMF is taken into account (Dundović & Sigl 2019); moreover, the substructures of the GMF themselves imprint smallscale anisotropies in cosmic rays (CRs) of petaelectronvolt (PeV) energies (Giacinti & Sigl 2012), and it is possible that such structures can be present at higher energies. Lastly, smallscale anisotropies can be detected separately in different regions of the sky, allowing us to probe, for example, different GMF structures independently.
Third, while most analyses have looked at the realspace correlation between UHECRs and the LSS (Kashti & Waxman 2008; Oikonomou et al. 2013; Abreu et al. 2010; Aab et al. 2015; Takami et al. 2009), we express our results here in terms of harmonicspace power spectra. These are common observables in cosmological studies and are based on a natural decomposition of the celestial sphere. They also allow for a straightforward visualisation of the main components of the astrophysical model (radial kernels, details of the galaxymatter connection), which is one of the main novel aspects of this work.
In this paper we introduce a formalism to model the AC and XC, and we apply it to a vanilla protononly model for UHECR injection in order to quantify the differences between the two observables and the detectability of the anisotropies on different scales with existing experimental facilities. We defer to upcoming work: (i) a more detailed discussion of the dependence of the XC on UHECR injection and source properties, (ii) a realistic treatment of the UHECR experimental setup, such as nonuniform sky coverage, and (iii) a full treatment of the effects of the GMF and xGMF on the signal.
This paper is organised as follows. In Sect. 2 we introduce the formalism to describe the UHECR flux, the distribution of galaxies, and the AC and XC. We apply this formalism to a hypothetical UHECR dataset in Sect. 3, where we obtain and compare the AC and the XC. We summarise our findings and conclude with an outlook for future work in Sect. 4. Appendix A presents useful standard formulae pertaining to angular power spectra.
2. Theoretical model
2.1. UHECR flux
Let ℰ(E_{inj}) be the (angleintegrated, isotropic) emissivity^{3} of CRs for a given galaxy (number of CRs of energy E_{inj} emitted per unit energy, per unit time):
The subscript ‘inj’ (injection) here indicates quantities evaluated in the rest frame of the emitting source. Due to the expansion of the Universe and to interactions between CRs and cosmic background light, the injected energy of a CR, whose energy at detection is E, is given by E_{inj}(E, z), where z is the redshift of the source. In the absence of scattering processes, the energy losses are adiabatic: E_{inj} = (1 + z)E. The differential emissivity (i.e., per unit solid angle) is ϵ := ℰ/4π, assuming isotropic emission. We parameterise the emissivity as a power law of energy:
Energies here are always expressed in exaelectronvolts (EeV) for convenience.
The quantity measured on Earth is the observed number of events per unit time, energy interval, detector area, solid angle on the sky, and (assuming source redshifts can be measured) redshift interval. We can relate this number to the emissivity through
where H(z) is the Hubble parameter and n_{s, c} is the volumetric number density of CR sources; we set c = 1, and we ignored subdominant lightcone and relativistic effects (Challinor & Lewis 2011; Bonvin & Durrer 2011).
We are interested in the number of UHECRs detected above a given energy threshold, E_{cut} (defined in the observer’s frame), and integrated over source redshifts, from the direction :
where χ(z) is the radial comoving distance.
We can write the number density of sources as , where δ_{s} is the galaxy overdensity. Assuming a nonevolving galaxy population, namely , and a powerlaw UHECR spectrum (as in Eq. (2)), we obtain:
2.1.1. Attenuation
The attenuation factor, α(E_{cut}, z; γ, Z), is defined as the number of events reaching the Earth with E > E_{cut}, divided by the number of events that would have reached the Earth if there were no energy losses at a given distance:
The attenuation, α, is a function of the energy cut and redshift as well as of the injection spectral slope and chemical composition. In terms of α, the directiondependent integral flux is
In this paper, to introduce the formalism, we chose to work with a toy protononly model with injection slope γ = 2.6 as in model (4) of d’Orfeuil et al. (2014) or model (i) of di Matteo & Tinyakov (2018). In order to obtain the attenuation factor for our injection model, we followed 10^{6} events with SimProp v2r4 (Aloisio et al. 2017) with energies above E = 10 EeV (with an upper cutoff of E = 10^{5} EeV), for redshifts up to z = 0.3, and counted the number of events reaching the Earth with E > E_{cut} for different values of E_{cut}. With SimProp we accounted for all energy losses, namely adiabatic losses and losses due to interactions with cosmic microwave background (CMB) photons and with extraGalactic background photons according to the model in Stecker et al. (2006). The UHECR radial kernels, defined in the next section, obtained from the attenuation factor, α, for different energies are shown in Fig. 1.
Fig. 1. Radial kernels for the two observables under consideration. The solid black line shows the approximate redshift distribution of galaxies in the 2MRS sample using the fit found by Ando et al. (2018). The red, yellow, and blue lines show the radial kernel for the UHECR flux (Eq. (10)) for the three energy thresholds studied here (40 EeV, 63 EeV, and 100 EeV, respectively). 
It should be noted that we assumed CR energy losses to be, to first order, isotropic, that is, we ignored angular anisotropies in the CMB and extraGalactic background light, which are completely negligible for our analysis. Moreover, for simplicity here we work with fullsky uniform coverage, but the analysis can be readily generalised to nonuniform and partial sky coverage.
2.1.2. Anisotropies
We defined the anisotropy in the UHECR distribution as the overdensity of rays detected as a function of sky position as
where is the skyaveraged UHECR flux. From the results in the previous section, this quantity is related to the threedimensional overdensity of UHECR sources through
where the UHECR radial kernel (or window function) is
Figure 1 shows the radial kernels for UHECRs with energy thresholds E_{cut} = 40, 63, and 100 EeV; as expected, the lower the energy, the farther UHECRs propagate.
2.2. Galaxies
We considered the AC of the UHECR anisotropy, Eq. (9), and its XC with the galaxy number count fluctuations. In particular, we worked with the projected overdensity of sources for a given galaxy sample,
where is the number of galaxies in a given direction () and is its average over the celestial sphere. This is related to the threedimensional galaxy overdensity via
where ϕ_{g}(χ) is the weighted distribution of galaxy distances. In general, we have assumed that we have redshift information for all galaxies in the catalogue and that we can use that information to apply a distancedependent weight, w(χ). In that case, the galaxy overdensity kernel, ϕ_{g}(χ), is given by
where is the comoving number density of galaxies in the sample.
If no weights are applied – namely, w(χ) = 1 – then
where is the angular number density of galaxies (i.e., the number of galaxies per steradian).
Figure 1 shows the radial kernel for a lowredshift galaxy survey, modelled after the 2MASS Redshift Survey (2MRS; Huchra et al. 2012). This constitutes one of the most complete fullsky spectroscopic lowredshift surveys, and we use it as our fiducial galaxy sample in this paper. In this work we consider fullsky datasets for simplicity, but the generalisation of our results for an incomplete sky coverage is straightforward. In the case of a realistic setup based on 2MRS, a sky coverage around 70% will only degrade the signal by a factor of .
2.3. Power spectra
We are interested in detecting the intrinsic anisotropies in the distribution of UHECRs by considering the different twopoint functions built from Δ_{CR} and Δ_{g}. A given observation of any of these fields will consist of both signal S and noise N: (where a, b ∈ {CR, g}). Assuming signal and noise to be uncorrelated, the corresponding power spectra can be split into both components, namely
where 𝒮_{ℓ} and 𝒩_{ℓ} are the power spectra of S and N, respectively. In our case, the signal is the intrinsic clustering of both UHECRs and galaxies due to the underlying LSS, while the noise is sourced by the discrete nature of both tracers as Poisson noise. A brief review of the mathematics behind angular power spectra is given in Appendix A.
2.3.1. Signal power spectra
The angular power spectrum, , between two projected quantities, Δ_{a} and Δ_{b}, is related to their threedimensional power spectrum, P_{ab}(z, k), by
where ϕ_{a} and ϕ_{b} are the radial kernels of both quantities.
The final piece of information needed in order to estimate the expected AC and XC signals is the power spectrum of the threedimensional overdensities δ_{s} and δ_{g}. In general, the clustering properties of galaxies and UHECR sources will depend on the specifics of the relationship between galaxies and dark matter, as well as on the astrophysical properties of the UHECR sources. To simplify the discussion, here we have assumed that all UHECR sources are also galaxies of the 2MASS sample (i.e., δ_{s} = δ_{g}).
At this point, one might be tempted to use a linear bias prescription (Mo & White 1996) to relate the galaxy and matter power spectra. However, as we show in Sect. 3, since the UHECR radial kernel peaks at z = 0 and covers only low redshifts, the CR flux AC probes mostly subhalo scales for which a nonperturbative description of structure formation is necessary. To achieve this, we used here a halo model prescription (Peacock & Smith 2000) that is based on the halo occupation distribution model used by Ando et al. (2018) to describe the 2MRS sample. In this model, the galaxy power spectrum is given by two contributions,
which are the socalled onehalo and twohalo terms. The former dominates on small scales and describes the distribution of galaxies within the halo, while the latter is governed by the clustering properties of dark matter haloes. The halo occupation distribution is then based on a prescription to assign central and satellite galaxies to haloes of different masses. Although we summarise this model in Appendix C, we refer the reader to Ando et al. (2018) and references therein for further details about the specifics of the halo occupation distribution model used.
2.3.2. Shot noise
Both projected overdensities, Δ_{CR} and Δ_{g}, are associated with discrete point processes, represented by the angular positions of the UHECRs and the galaxies in each sample. In this case, even in the absence of intrinsic correlations between the different fields, their power spectra receive a nonzero white contribution, given by
where () is the angular number density of points in sample a or b, and is the angular number density of points shared in common. In our case this corresponds to the number of UHECRs originating from galaxies in the galaxy sample. For simplicity, we have assumed that the galaxy survey under consideration is sufficiently complete, such that all UHECRs are associated with an observed galaxy. In this case, the shot noise contributions to the power spectra are
Since typically , then , and therefore we neglect in what follows. We have explicitly checked that the crossnoise can indeed be safely neglected in all our estimates and numerical results.
We note that when nonflat weights are applied to the galaxy catalogue, the resulting noise power spectrum reads
For w(χ) = 1, Eq. (14) holds, and we recover the result in Eq. (19).
2.3.3. Optimal weights
We can use the results in this section to derive optimal weights w(χ) to maximise the signaltonoise of the galaxyUHECR XC. We pixelised the celestial sphere and considered the UHECRs in a given pixel p, Φ_{p}, as well as the vector N_{p, i}, which contains the number of galaxies along the same pixel in intervals of distance χ_{i}. The optimal weights w_{i} := w(χ_{i}) can be found by maximising the likelihood of Φ_{p} given N_{p, i} (Alonso et al. 2020) and are given by the socalled Wiener filter, namely
Here, Cov(x, y) is the covariance matrix of two vectors, x an y.
Assuming Poisson statistics, we can use the results from the previous section to show (see Appendix B) that
In hindsight, this result is obvious: By inspecting Eqs. (10) and (13), we see that the optimal weights modify the radial galaxy kernel, ϕ_{g}, to make it identical to the UHECR kernel, ϕ_{CR}, thereby building the most likely estimate of the UHECR flux map from the galaxy positions. As we will see, this involves upweighting galaxies at low redshifts, where the UHECRs that reach the Earth most likely originate, but few galaxies can be found due to volume effects. We notice that the weights are completely driven by the theoretical model for UHECR propagation and do not depend on the actual data (or their errors). We show how the use of optimal weights can improve the S/N for the XC in Sect. 3.2.
2.4. Intervening magnetic fields
The Milky Way is host to a magnetic field of a few μG (Boulanger et al. 2018), which is the screen that befogs UHECR sources. The variety of the parametric models of the GMF – which disagree on the GMF functional forms and parameters, particularly on GMF substructures – reflects the complexity of the GMF, and, at the moment, none can be taken at face value (Boulanger et al. 2018; Unger & Farrar 2017). As a guideline, we can expect the GMF to deflect a UHECR with energy E = 100 EeV by a few degrees for most parts of the sky, except for certain directions close to the Galactic plane (d’Orfeuil et al. 2014; see also Pshirkov et al. 2013). UHECRs will also be affected by any intervening xGMF, whose strength, shape, and filling factors vary by several orders of magnitude for different models and estimates (Subramanian 2016); however, the xGMF is believed to have a subdominant effect on largescale UHECR propagation (Pshirkov et al. 2016).
Simple prescriptions to account for part of the effects related to the GMF and the xGMF include smearing the map of sources below a certain angular scale or mixing that map with an isotropic one (similar to what is done to take into account catalogue incompleteness beyond a certain distance; see e.g., Koers & Tinyakov 2009). Smearing the source map in our language is as simple as introducing a (Gaussian) beam in the signals as
for the AC and XC, respectively, where
and θ_{smear} is the smearing angle.
However, these solutions tend to be rather artificial and inaccurately destroy potential signals or structures in the spectra that we are examining. More precisely, we know that the largest effects due to the smallscale GMF are not isotropic and in fact vary quite considerably across the sky (see e.g., Pshirkov et al. 2013; di Matteo & Tinyakov 2018). More precisely, Pshirkov et al. (2013) found that the UHECR deflections are majorated by the function
where b is elevation. Therefore, if we smear the whole sky with the same smearing angle, we are not faithfully representing the sky, and, depending on the smearing angle, we either underestimate the deflections in certain regions or overestimate them in other regions (or both). Moreover, these solutions do not account for the largescale galactic field, which is the dominant effect and cannot be described by a simple smearing.
For these reasons, and in order to best introduce the method, in this first theoretical work we take a pragmatic approach and, keeping in mind all the caveats listed above, only briefly discuss the effect of a (constant) smearing angle on the AC and XC (see also Appendix D). We neglect all other effects of intervening magnetic fields.
3. Results
3.1. Signaltonoise ratio
We estimate the S/N of the UHECR anisotropies as the square root of the Fisher matrix element corresponding to an effective amplitude parameter, A_{CR}, multiplying the signal component of Δ_{CR} with a fiducial value A_{CR} = 1 (Heavens 2009), namely
where S_{ℓ} is a vector containing the signal contribution to the power spectra under consideration, Cov is the covariance matrix of those power spectra, and S/N_{ℓ} is the S/N of a single ℓ mode. If the fields being correlated are Gaussian (Δ_{CR}, Δ_{g} in our case), the covariance matrix can be estimated using Wick’s theorem to be
where Δℓ is the size of the multipole bin.
At this point we consider three different cases. The first is AC only. In this case we only have a measurement of the UHECR AC, . The S/N is given by
The second is XC only. In this case we only use the XC, . The S/N is given by
In the third case we use all available data, that is, a data vector . Although this is the manifestly optimal scenario, XCs are arguably safer than ACs in terms of systematic errors, and therefore it is interesting to quantify the loss of information if only XCs are used.
Studying these three cases allows us to explore the benefits of using XCs versus ACs, as well as the relative amount of information in each of the different twopoint functions. Given the relatively small number of UHECRs currently measured, shot noise in the UHECR flux is the dominant contribution to the uncertainties. Comparing Eqs. (30) and (31), we can see that the S/N scales as and for cases 1 and 2, respectively, highlighting the potential of XCs for achieving a detection.
3.2. Power spectra and signaltonoise
The energy at which we chose to cut the UHECR integral spectrum, E_{cut}, determines the UHECR propagation horizon, which in turns determines the strength of the anisotropy. Moreover, the choice of E_{cut} for a given UHECR spectrum also determines the number of UHECR events we have to sample the anisotropic angular distribution. We expect a tradeoff between the two. At low energies the UHECR sample contains many more events than at high energies because the UHECR spectrum drops very steeply with energy. However, for the range of energies we are interested in, the galaxy sample is much larger; as such, this does not have as strong an effect for the XC as it does for the traditional AC (whose noise is determined by the number of UHECR events). Moreover, at low energies UHECRs propagate farther, and the larger lineofsight averaging can dilute the expected anisotropy. Lastly, the effects of intervening magnetic fields are stronger – this is expected to have a significant impact on the anisotropies, albeit less so for the XC compared to the AC thanks to its stability against systematics. At high energies the UHECR horizon is smaller, UHECRs undergo smaller deflections, and the anisotropy should be more pronounced, but the number of events drops dramatically.
In order to determine at which energy we have the best chances of detecting the XC, we chose to work with three energy cuts at: E_{cut} = 10^{19.6} eV ≃ 40 EeV, E_{cut} = 10^{19.8} eV ≃ 63 EeV, and E_{cut} = 10^{20} eV = 100 EeV. In a realistic scenario, based on currently available data (Alves Batista et al. 2019), we can expect to have about N_{CR} = 1000, N_{CR} = 200, and N_{CR} = 30 over the full sky for these three energy cuts, respectively. While this does not fully reflect a realistic situation – mostly because of magnetic deflections, which we do not take into account, and because current experimental facilities are limited in their field of view – our results nonetheless present a fair comparison between the two measures (AC and XC). This is owing to the fact that all the salient information regarding UHECR datasets is represented in our estimates, namely the energy cut and with it all the UHECR energy losses, the available or expected number of events at that energy, and the angular resolution representative of what current experiments can do.
In Fig. 2 we show the expected signal for the AC (left panel) and the XC (right panel). Colours refer to the three energy cuts discussed above, namely red for E_{cut} ≃ 40 EeV, yellow for E_{cut} ≃ 63 EeV, and blue for E_{cut} = 100 EeV. The dashed and dotted curves show the onehalo and twohalo contributions to the total power spectrum, with the sum of both shown by the solid curves. For simplicity, we have not included any beam smoothing in the plot. We can see how the signal for the XC is lower than for the AC, as is expected from the fact that the XC mixes two different radial kernels. If we employed optimal weights for the XC, the signal would become identical to that of the AC. In our simplistic linear treatment of perturbations, this happens because the UHECR and galaxy kernels would be identical. The statistical uncertainties for both correlation functions, however, would be different given their different shot noise levels.
Fig. 2. Angular AC and XC power spectra considered in this work. Dotted and dashed curves respectively refer to the one and twohalo contribution to the total signal (solid curves). 
To better understand the role of the different uncertainties on the theoretical signal, in Fig. 3 we again show the expected signal as in Fig. 2 (solid curves, same colour code) and include a 1° Gaussian smoothing beam to account for the angular resolution of UHECR experiments (for reference, we also show the beamfree prediction as dashed lines). On top of it, we present the corresponding ℓbinned 1σ error bars as shaded boxes for 20 logspace multipole bins between ℓ_{min} = 2 and ℓ_{max} = 1000. If we compare the leftmost and central panels, namely AC versus XC, it is easy to see how the range of multipoles where error bars are small enough to allow a detection is larger for XCs than for ACs for the E_{cut} ≃ 40 EeV and E_{cut} ≃ 63 EeV cases. However, for the sparser UHECR sample with E_{cut} = 100 EeV, the opposite applies; more precisely, the detectable range of multipoles for the XC is smaller and pushed towards higher ℓ compared to the AC. This is due to a combination of two factors: for the higher end of UHECR energies, the propagation horizon of UHECRs is small and the UHECR sky looks more anisotropic, boosting the AC. At the same time, the mismatch in kernels is prominent, the more so the higher the energy, and this drives the XC signal down. Combined with the larger shot noise in the UHECR data, this can explain the performance of the 100 EeV case – indeed, the UHECR shot noise is the main factor that prevents a detection of the signal at midℓ values (the perℓ signal is 1σ compatible with zero; see below).
Fig. 3. Expected power spectra and ℓbinned 1σ uncertainties (shaded boxes), including a 1° Gaussian smoothing beam to account for the angular resolution of UHECR experiments (solid curves). For reference, horizontal lines in the leftmost plots denote shot noise levels, and the dashed curves show the beamfree prediction. 
In the rightmost panel of Fig. 3 we show the XC signal when we apply theoretical optimal weights. In this case the highestenergy set performs the best, and this is expected from the previous arguments: the signal is boosted back up to the same level of the AC because the kernels of galaxies and UHECRs now coincide. Additionally, while the uncertainty increases with energy as both samples become sparser, it is not large enough to hide the XC signal. It is worth noticing that the increase in galaxy power that we expect towards lower redshifts is significantly less relevant than the matching of the radial kernels.
In practice, using optimal weights may not be possible given the uncertainties in the radial kernel for UHECRs (we do not yet know the actual injection spectrum). The availability of redshift information in the galaxy catalogue, however, would allow us to turn this into an advantage: the UHECR kernel could be reconstructed by modifying the galaxy weights to maximise the signaltonoise, essentially following the ‘clustering redshifts’ method used to reconstruct unknown redshift distributions in weak lensing data (Newman 2008).
To quantify the improvement in detectability brought by the XC, in Fig. 4 we present the cumulative S/N for all the data combinations discussed in Sect. 3.1, viz. AC alone (leftmost panel), XC alone (central panel), and all the data combined into a single data vector S_{ℓ} (rightmost panel). In each panel, the left half shows the cumulative S/N as a function of the maximum multipole, ℓ_{max}, whilst the right half is for the cumulative S/N as a function of the minimum multipole, ℓ_{min}. In both cases, the case with all the data combined has, unsurprisingly, the largest S/N, but the contributions from AC and XC come from different angular scales. These scales in turn are sensitive to different redshift ranges depending on E_{cut}, which sets the propagation depth for UHECRs. This highlights the complementarity of the two observables.
Fig. 4. S/N for UHECR flux anisotropies from different combinations of data, namely UHECR AC in the leftmost panel, XC in the central panel, and the combination of all data in the rightmost panel. In each panel, the left half shows the cumulative S/N as a function of the maximum multipole, ℓ_{max}, whereas the right half is for the cumulative S/N as a function of the minimum multipole, ℓ_{min}. The horizontal dashed line marks the 3σ threshold for detection. 
The aforementioned sensitivity to different angular scales can be captured better by looking at Fig. 5, where we show the contribution to the total S/N from each integer multipole, S/N_{ℓ}. The colour code is the same as throughout the paper, and we mark with horizontal dashed lines the thresholds corresponding to 1, 2, and 3σ evidence for a oneparameter amplitude fit. These panels can be interpreted as evidence for anisotropy on a given scale; it is clear that the XC with galaxies helps to push the detectability of the signal to smaller scales (i.e. larger ℓ values). This perℓ S/N_{ℓ} is a useful quantity for assessing whether the AC or the XC is the best observable for detecting the anisotropy in UHECRs, assuming that UHECRs trace the LSS.
Fig. 5. S/N per multipole, S/N_{ℓ}, for the AC signal, the XC signal with both normal and optimal weights, and their combination, AC+XC (leftmost, central, and rightmost panel, respectively). Different colours refer to different energy cuts, and the three horizontal dashed lines show the thresholds for 1, 2, and 3σ detection. 
The sensitivity of the XC to smallscale anisotropies can be precious in realistic situations for two further reasons. First, a single Earthbased experiment is blind to a large fraction of the sky (roughly speaking, one celestial hemisphere); galaxy catalogues can also have incomplete sky coverage. Moreover, it might be advantageous (see our discussion of the directiondependent magnetic deflections in Sect. 2.4) to restrict the UHECR dataset to a portion of the sky to maximise the chances for a clean detection. In all these situations, the low harmonic multipoles are the most affected by these sky cuts. Second, if two experiments join their datasets, as the Telescope Array and Pierre Auger collaborations have done in their harmonic AC analysis, they need to crosscalibrate their sets, and this crosscalibration introduces errors that are significantly larger for low multipoles than for high multipoles (AbuZayyad et al. 2013; Aab et al. 2014; Deligny 2015; di Matteo et al. 2020).
As we argued in Sect. 2.4, there is no shortcut to account for the effects of the GMF on the AC and XC. Nonetheless, it is instructive to look at how the signal degrades with a simple smearing of the source map. To this end, we have plotted in Fig. 6 the total S/N for ℓ = [2, 1000] as a function of the smearing angle, θ_{smear}, of the galaxy map for the same energy cuts we have used so far. According to Eq. (26), the deflections for 40 EV rigidity peak at around 7° near the Galactic centre, whereas more than half of the sky would be well described by a 2.5° smearing. It should be noted that, as mentioned in the introduction, the composition of UHECRs at the highest energies is unknown (Castellina 2020; Bergman 2019), and a heavier composition would imply lower rigidities and larger deflections. The smearing impacts the highmultipole regions in the XC more than it does in the AC, as expected, and degrades the XC more prominently at larger smearing angles (see also Appendix D).
Fig. 6. Total S/N, , as a function of the smearing angle for the AC signal, the XC signal with both normal and optimal weights, and their combination, AC+XC (leftmost, central, and rightmost panel, respectively). Different colours refer to different energy cuts, and the horizontal dashed line shows the thresholds for 3σ detection. 
The XC is strongly dependent on energy and the choice of weights. This means that we can disentangle smallscale anisotropies caused by the propagation through the GMF from the intrinsic anisotropies inherited from the LSS. In particular, any GMFinduced signal is erased at higher energies because UHECRs go straighter, whereas the LSSinherited signal is enhanced because of the smaller propagation horizon. Moreover, any GMFinduced signal is indifferent to the weights we apply, whereas the signal from LSS anisotropies is strongly enhanced with the use of optimal weights.
Before closing this section, we remark that in a real experiment there will be modelled and unmodelled systematic errors to take into account. Systematic errors are expected to contribute to the AC more significantly than to the XC, particularly on large scales (lowℓ end), for example the crosscalibration of two UHECR datasets. On the other hand, biassed redshift information in galaxy catalogues or UHECR injection properties will affect both observables. To be clear, while the galaxy catalogue and the optimal weights, which depend on UHECR data for the reconstruction of the injection properties, do not enter into the AC obtained from UHECR data alone, they are needed once we test the source model (e.g., so that UHECRs correlate with the LSS). Hence, once systematic effects are taken into account, the S/N for the AC may decrease more than the S/N for the XC. This is one further motivation to explore the possibilities of and improvements from the use of XCs in UHECR anisotropy studies.
4. Conclusions and outlook
In this work we have introduced a new observable for UHECR physics: the harmonicspace XC between the arrival directions of UHECRs and the distribution of the cosmic LSS as mapped by galaxies (Eq. (16)). We have developed the main theoretical tools necessary to model the signal and its uncertainties.
The takeaway points of this study are as follows.

The XC can be easier to detect than the UHECR AC for a range of energies and multipoles (see Figs. 3 and 5). This performance is mostly driven by the sheer number of galaxies that can trace the underlying LSS distribution, which is assumed to be the baseline distribution for both the UHECR flux and the galaxy angular distribution.

The XC is more sensitive to smallscale angular anisotropies than the AC; conversely, the AC is more sensitive to largescale anisotropies. This finding could therefore be instrumental in understanding properties of UHECR sources that would not be accessible otherwise.

It is in principle possible to optimise the XC signal by assigning optimal redshiftdependent weights to sources in the galaxy catalogue in order to match the UHECR radial kernel as determined by UHECR energy losses. Since matching the kernels has a strong impact on the XC, it could be possible to use this effect to reverse engineer the injection model (which defines the radial kernel).

The great disruptor of UHECR anisotropies is the GMF. The XC, with its higher S/N and sensitivity to small angular scales, could be very useful in understanding the properties of the GMF (although we have not explored this angle here). Moreover, it may be possible, in the near future, to exploit a tomographic approach to disentangle the effects of intervening magnetic fields from different injection spectra and study different regions of the sky separately.
In our treatment we do not take any experimental uncertainties into account, with the exception of the experimental UHECR angular resolution. Moreover, we limit ourselves to a protononly injection model and do not include the effects of the intervening magnetic fields. This choice was made in order to underline the physics behind our proposal and method. This method can be readily generalised and extended to include the (theoretical and experimental) properties of the different galaxy and UHECR catalogues and the different injection models as well as to separate the number of events and energy cuts in order to best forecast the possibilities of present and upcoming UHECR datasets.
Moreover, in this initial work we have made the case for the XC between UHECRs and galaxies, but the logic and methods we have developed can be applied to other XCs with different matter tracers and different messengers. The distribution of visible matter in the sky can be traced not only by galaxies, but also by the thermal SunyaevZeldovich effect. This effect is produced by the inverse Compton scattering of CMB photons by hot electrons along the line of sight. Because a thermal SunyaevZeldovich map is a map of CMB photons, it is very accurate down to angles much smaller than a degree, and its signal peaks at low redshifts (Erler et al. 2018). This XC could therefore be useful in further disentangling the astrophysical properties of UHECR sources.
Charged UHECRs are not the only highenergy messengers whose production mechanisms and sources are unknown. The IceCube collaboration has recently detected a few highenergy astrophysical neutrinos with energies above 1 PeV (Aartsen et al. 2014). Such neutrinos are expected to be produced in the same extreme astrophysical sources as UHECRs and/or in their immediate surroundings. The XCs between neutrinos and the LSS will then inform us about the properties of the highestenergy astrophysical engines (see Fang et al. 2020 as well as Ando et al. 2015). Without the use of XCs (due to the very small number of neutrino events in current data), the detection of the anisotropic pattern could remain challenging into the foreseeable future (Fang et al. 2020; Sapienza 2020; Allison et al. 2016; Nelles 2019). Since neutrinos interact extremely weakly, they can propagate unhampered for long distances: their horizon is almost the entire visible Universe. Therefore, complementary information could be extracted from crosscorrelating neutrinos with other tracers (in addition to galaxies), including CMB lensing (Planck Collaboration VIII 2020) and cosmic shear surveys (Mandelbaum 2018), both of which trace the overall matter distribution in the Universe, including both its dark and luminous components, out to higher redshifts with broader kernels (see Fornengo et al. 2015; Cuoco et al. 2015; Branchini et al. 2017; Ammazzalorso et al. 2020 for the analogous analysis with γ rays). Measuring these XCs could reveal whether the most energetic particle accelerators in the Universe preferentially reside in highdensity visible or dark environments.
Acknowledgments
FU wishes to thank A. di Matteo for useful correspondence, and P. Tinyakov for valuable comments on the manuscript. FU is supported by the European Regional Development Fund (ESIF/ERDF) and the Czech Ministry of Education, Youth and Sports (MEYS) through Project CoGraDS – CZ.02.1.01/0.0/0.0/15_003/0000437. SC is supported by the Italian Ministry of Education, University and Research (MIUR) through Rita Levi Montalcini project ‘PROMETHEUS – Probing and Relating Observables with Multiwavelength Experiments To Help Enlightening the Universe’s Structure’, and by the ‘Departments of Excellence 20182022’ Grant awarded by MIUR (L. 232/2016). DA acknowledges support from the Beecroft Trust, and from the Science and Technology Facilities Council through an Ernest Rutherford Fellowship, grant reference ST/P004474/1.
References
 Aab, A., Abreu, P., Aglietta, M., et al. 2014, ApJ, 794, 172 [CrossRef] [Google Scholar]
 Aab, A., Abreu, P., Aglietta, M., et al. 2015, ApJ, 804, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Aab, A., Abreu, P., Aglietta, M., et al. 2017, Science, 357, 1266 [NASA ADS] [CrossRef] [Google Scholar]
 Aab, A., Abreu, P., Aglietta, M., et al. 2019, ICRC, C4, 20 [Google Scholar]
 Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, Phys. Rev. Lett., 113 [CrossRef] [Google Scholar]
 Abbasi, R. U., AbuZayyad, T., Allen, M., et al. 2008, Phys. Rev. Lett., 100 [CrossRef] [Google Scholar]
 Abbasi, R. U., AbuZayyad, T., Allen, M., et al. 2014, ApJ, 790, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Abbasi, R. U., AbuZayyad, T., Allen, M., et al. 2017, Astropart. Phys., 86, 21 [CrossRef] [Google Scholar]
 Abraham, J., Abreu, P., Aglietta, M., et al. 2008, Phys. Rev. Lett., 101 [CrossRef] [Google Scholar]
 Abreu, P., Aglietta, M., Ahn, E. J., et al. 2010, Astropart. Phys., 34, 314 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Abreu, P., Aglietta, M., Ahn, E. J., et al. 2012, Astropart. Phys., 35, 354 [CrossRef] [Google Scholar]
 AbuZayyad, T., Allen, M., Anderson, R., et al. 2013, 33rd International Cosmic Ray Conference [Google Scholar]
 Allison, P., Bard, R., Beatty, J. J., et al. 2016, Phys. Rev. D, 93 [CrossRef] [Google Scholar]
 Aloisio, R., Boncioli, D., Di Matteo, A., et al. 2017, JCAP, 1711, 009 [CrossRef] [Google Scholar]
 Alonso, D., Salvador, A. I., Sánchez, F. J., et al. 2015, MNRAS, 449, 670 [NASA ADS] [CrossRef] [Google Scholar]
 Alonso, D., Cusin, G., Ferreira, P. G., & Pitrou, C. 2020, Phys. Rev. D, 102, 023002 [CrossRef] [Google Scholar]
 Alvarez, E., Cuoco, A., Mirabal, N., & Zaharijas, G. 2016, JCAP, 12, 023 [CrossRef] [Google Scholar]
 Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, Front. Astron. Space Sci., 6, 23 [CrossRef] [Google Scholar]
 Ammazzalorso, S., Gruen, D., Regis, M., et al. 2020, Phys. Rev. Lett., 124 [CrossRef] [Google Scholar]
 Ando, S., Tamborra, I., & Zandanel, F. 2015, Phys. Rev. Lett., 115 [CrossRef] [Google Scholar]
 Ando, S., BenoitLévy, A., & Komatsu, E. 2018, MNRAS, 473, 4318 [CrossRef] [Google Scholar]
 Bergman, D. 2019, PoS, ICRC2019, 190 [Google Scholar]
 Bonvin, C., & Durrer, R. 2011, Phys. Rev. D, 84 [Google Scholar]
 Boulanger, F., Enßlin, T., Fletcher, A., et al. 2018, JCAP, 1808, 049 [Google Scholar]
 Branchini, E., Camera, S., Cuoco, A., et al. 2017, ApJS, 228, 8 [Google Scholar]
 Caccianiga, L. 2020, PoS, ICRC2019, 206 [Google Scholar]
 Camera, S., Fornasa, M., Fornengo, N., & Regis, M. 2013, ApJ, 771, L5 [Google Scholar]
 Castellina, A. 2020, PoS, ICRC2019, 004 [Google Scholar]
 Challinor, A., & Lewis, A. 2011, Phys. Rev. D, 84 [Google Scholar]
 Cooray, A., & Sheth, R. K. 2002, Phys. Rept., 372, 1 [Google Scholar]
 Cuoco, A., Xia, J.Q., Regis, M., et al. 2015, ApJS, 221, 29 [Google Scholar]
 Deligny, O. 2015, PoS, ICRC2015, 395 [Google Scholar]
 di Matteo, A., & Tinyakov, P. 2018, MNRAS, 476, 715 [Google Scholar]
 di Matteo, A., Deligny, O., Kawata, K., et al. 2018, JPS Conf. Proc., 19 [Google Scholar]
 di Matteo, A., Aloisio, R., Boncioli, D., et al. 2020, PoS, ICRC2019, 439 [Google Scholar]
 d’Orfeuil, B. R., Allard, D., Lachaud, C., et al. 2014, A&A, 567, A81 [CrossRef] [EDP Sciences] [Google Scholar]
 Dundović, A., & Sigl, G. 2019, JCAP, 1901, 018 [Google Scholar]
 Erler, J., Basu, K., Chluba, J., & Bertoldi, F. 2018, MNRAS, 476, 3360 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, K., Banerjee, A., Charles, E., & Omori, Y. 2020, ApJ, 894, 112 [Google Scholar]
 Fornengo, N., & Regis, M. 2014, Front. Phys., 2, 6 [Google Scholar]
 Fornengo, N., Perotto, L., Regis, M., & Camera, S. 2015, ApJ, 802, L1 [Google Scholar]
 Giacinti, G., & Sigl, G. 2012, Phys. Rev. Lett., 109 [Google Scholar]
 Greisen, K. 1966, Phys. Rev. Lett., 16, 748 [Google Scholar]
 Heavens, A. 2009, Statistical Techniques in Cosmology [Google Scholar]
 Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, ApJS, 199, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Kashti, T., & Waxman, E. 2008, JCAP, 05, 006 [Google Scholar]
 Koers, H. B. J., & Tinyakov, P. 2009, JCAP, 0904, 003 [Google Scholar]
 Kotera, K., & Olinto, A. V. 2011, ARA&A, 49, 119 [Google Scholar]
 Linsley, J., Scarsi, L., & Rossi, B. 1961, Phys. Rev. Lett., 6, 485 [Google Scholar]
 Mandelbaum, R. 2018, ARA&A, 56, 393 [Google Scholar]
 Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 [Google Scholar]
 Nelles, A. 2019, EPJ Web Conf., 216, 01008 [Google Scholar]
 Newman, J. A. 2008, ApJ, 684, 88 [Google Scholar]
 Oikonomou, F., Connolly, A., Abdalla, F. B., et al. 2013, JCAP, 05, 015 [Google Scholar]
 Pan, J., & Coles, P. 2000, MNRAS, 318, L51 [Google Scholar]
 Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144 [NASA ADS] [CrossRef] [Google Scholar]
 Pinetti, E., Camera, S., Fornengo, N., & Regis, M. 2020, JCAP, 07, 044 [Google Scholar]
 Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2013, MNRAS, 436, 2326 [Google Scholar]
 Planck Collaboration VIII. 2020, A&A, 641, A8 [CrossRef] [EDP Sciences] [Google Scholar]
 Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, Phys. Rev. Lett., 116 [Google Scholar]
 Sapienza, P. 2020, J. Phys. Conf. Ser., 1342 [Google Scholar]
 Scrimgeour, M., Davis, T., Blake, C., et al. 2012, MNRAS, 425, 116 [Google Scholar]
 Sommers, P. 2001, Astropart. Phys., 14, 271 [Google Scholar]
 Stecker, F. W., Malkan, M. A., & Scully, S. T. 2006, ApJ, 648, 774 [Google Scholar]
 Subramanian, K. 2016, Rept. Prog. Phys., 79 [Google Scholar]
 Takami, H., Nishimichi, T., Yahata, K., & Sato, K. 2009, JCAP, 06, 031 [Google Scholar]
 Tinyakov, P. G., & Urban, F. R. 2015, J. Exp. Theor. Phys., 120, 533 [Google Scholar]
 Tinyakov, P. G., Urban, F. R., Ivanov, D., Thomson, G. B., & Tirone, A. H. 2016, MNRAS, 460, 3479 [Google Scholar]
 Unger, M., & Farrar, G. R. 2017, 35th Int. Cosmic Ray Conf., 35, 558 [Google Scholar]
 Zatsepin, G. T., & Kuzmin, V. A. 1966, JETP Lett., 4, 78 [Pisma Zh. Eksp. Teor. Fiz. 4,114(1966)] [NASA ADS] [Google Scholar]
Appendix A: Power spectra
Threedimensional fields δ_{a}(x) can be decomposed into their Fourier modes
whose covariance is the power spectrum P_{ab}(k). Assuming statistical homogeneity and isotropy, P_{ab}(k) is implicitly defined by
where the angle brackets denote averaging over ensemble realisations of the random fields inside them.
Equivalently, twodimensional fields can be decomposed into their harmonic coefficients:
where Ω = (θ, φ) is the solid angle on the sky, Y_{ℓm} are the spherical harmonic functions, is the lineofsight direction, and in this work a = {CR, g}. The covariance of the Δ_{ℓm} is the angular power spectrum, , defined as
For two projected fields, Δ_{a} and Δ_{b}, associated with threedimensional fields δ_{a} and δ_{b} via radial kernels ϕ_{a} and ϕ_{b} (as in Eqs. 9 and 12), their Fourier and harmonicspace power spectra are related through
where j_{ℓ} is the spherical Bessel function of order ℓ. For broad kernels, we can use the Limber approximation, , in which case the previous relation simplifies to Eq. (16).
Appendix B: Optimal weights
Here we derive the choice of optimal weights discussed in Sect. 2.3.3. The derivation is a standard result in statistics and follows the discussion in Appendix A of Alonso et al. (2020).
We consider a vector of N measurements x = (x_{1}, …, x_{N}) and the problem of finding the linear combination of this vector that provides the best estimator of a given quantity, y. Assuming Gaussian statistics, the conditional probability is
where z ≡ (y, x_{1}, …, x_{N}), and C_{ab} is the covariance matrix between a and b^{4}. The minimumvariance estimator for y given this distribution coincides with its mean, which is given by
The linear coefficients w are the socalled Wiener filter.
If y is the UHECR flux and x is a set of galaxy overdensity maps at different radial shells with comoving width dχ, in the shot noisedominated regime C_{xy} and C_{xx} are given by
where we have ignored all rindependent prefactors. Therefore, the Wiener filter in this case is
Appendix C: Halo occupation distribution
Halo occupation distribution models have been used profusely in the literature to describe the connection between the galaxy number density and the matter overdensity. Here we give a brief overview of the specifics of the model used in this work to describe the lowredshift 2MRS sample, which follows Ando et al. (2018).
Within the halo model (Peacock & Smith 2000; Cooray & Sheth 2002), all matter in the Universe can be found in haloes of different masses, and therefore the fluctuations of a given quantity, x, can be described in terms of its distribution around haloes as a function of halo mass, u_{x}(r, M) (also called the halo profile), and the correlated distribution of haloes on large scales. In this formalism, the power spectrum between two quantities, x and y, receives two contributions, coming from the correlations between mass elements belonging to the same halo and mass elements in different haloes (the socalled onehalo and twohalo terms), as , with
where n_{h}(M) and b_{h}(M) are respectively the halo mass function and the halo bias, P_{lin}(k) is the linear matter power spectrum, and u_{x}(k, M) is the Fourier transform of the halo profile.
In our case, we want to model the galaxy overdensity, δ_{g}, and therefore we need to specify the halo galaxy density profile. For this, we employed the formalism used in Ando et al. (2018):
where N_{c}(M) and N_{s}(M) are the number of central and satellite galaxies, the latter of which are distributed according to u_{s}(r). Centrals and satellites are distributed according to Bernoulli and Poisson distributions, respectively. Their mean values and the satellite profile are modelled as a function of mass as:
where Θ is the Heaviside function. The free parameters of the model are M_{min}, M_{1}, σ_{ log M}, r_{max}, g/r_{s}, r_{s, g}/r_{s}, and α, where r_{s} is the massdependent halo scale radius. We used the bestfit values found by Ando et al. (2018) for these parameters in our calculation.
Appendix D: Magnetic deflections
The effects of the GMF deflections are not the same across harmonic multipoles. We expect that, for both the AC and the XC, small scales would be more affected by the deflections. We illustrate this in Fig. D.1, where we show the equivalent of Fig. 6 but for multipoles in the four half decades: ℓ ∈ [3, 10[, ℓ ∈ [10, 33[, ℓ ∈ [33, 100[, and ℓ ∈ [100, 333[. As anticipated, the smearing suppresses the power at small scales more incisively, for both the AC and the XC, with the XC more affected. Nonetheless, in regions of the sky where the GMF is small, for example around the Galactic polar caps, the XC has better chances to be detected than the AC at large ℓ (small scales). In reading this figure one should keep at least three simplifications in mind: The larger but physically different effects of the largescale GMF are not included; the smearing angle is constant across the sky (see Eq. (26)); and the choice of a Gaussian smearing is arbitrary as other types of beams might reproduce the actual deflections more faithfully.
Fig. D.1. Total S/N per halfdecade, , with (ℓ_{min}, ℓ_{max}) in [3, 10[ (top row), [10, 33[ (second row), [33, 100[ (third row), and [100, 333[ (bottom row), as a function of the smearing angle for the AC signal, the XC signal with both normal and optimal weights, and their combination, AC+XC (leftmost, central, and rightmost panel, respectively, in each row). Different colours refer to different energy cuts, and the horizontal dashed line shows the thresholds for 3σ detection. 
All Figures
Fig. 1. Radial kernels for the two observables under consideration. The solid black line shows the approximate redshift distribution of galaxies in the 2MRS sample using the fit found by Ando et al. (2018). The red, yellow, and blue lines show the radial kernel for the UHECR flux (Eq. (10)) for the three energy thresholds studied here (40 EeV, 63 EeV, and 100 EeV, respectively). 

In the text 
Fig. 2. Angular AC and XC power spectra considered in this work. Dotted and dashed curves respectively refer to the one and twohalo contribution to the total signal (solid curves). 

In the text 
Fig. 3. Expected power spectra and ℓbinned 1σ uncertainties (shaded boxes), including a 1° Gaussian smoothing beam to account for the angular resolution of UHECR experiments (solid curves). For reference, horizontal lines in the leftmost plots denote shot noise levels, and the dashed curves show the beamfree prediction. 

In the text 
Fig. 4. S/N for UHECR flux anisotropies from different combinations of data, namely UHECR AC in the leftmost panel, XC in the central panel, and the combination of all data in the rightmost panel. In each panel, the left half shows the cumulative S/N as a function of the maximum multipole, ℓ_{max}, whereas the right half is for the cumulative S/N as a function of the minimum multipole, ℓ_{min}. The horizontal dashed line marks the 3σ threshold for detection. 

In the text 
Fig. 5. S/N per multipole, S/N_{ℓ}, for the AC signal, the XC signal with both normal and optimal weights, and their combination, AC+XC (leftmost, central, and rightmost panel, respectively). Different colours refer to different energy cuts, and the three horizontal dashed lines show the thresholds for 1, 2, and 3σ detection. 

In the text 
Fig. 6. Total S/N, , as a function of the smearing angle for the AC signal, the XC signal with both normal and optimal weights, and their combination, AC+XC (leftmost, central, and rightmost panel, respectively). Different colours refer to different energy cuts, and the horizontal dashed line shows the thresholds for 3σ detection. 

In the text 
Fig. D.1. Total S/N per halfdecade, , with (ℓ_{min}, ℓ_{max}) in [3, 10[ (top row), [10, 33[ (second row), [33, 100[ (third row), and [100, 333[ (bottom row), as a function of the smearing angle for the AC signal, the XC signal with both normal and optimal weights, and their combination, AC+XC (leftmost, central, and rightmost panel, respectively, in each row). Different colours refer to different energy cuts, and the horizontal dashed line shows the thresholds for 3σ detection. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.