A GLOBAL AUTOCORRELATION STUDY AFTER THE FIRST AUGER DATA: IMPACT ON THE NUMBER DENSITY OF UHECR SOURCES

Evidence is now emerging that ultrahigh energy cosmic rays (UHECRs) have an astrophysical origin, as opposed to being generated in exotic top-down models: the detection of a spectral suppression consistent with the Greisen–Zatsepin–Kuzmin (GZK) effect (Greisen 1966; Zatsepin & Kuzmin 1966) by both HiRes (Abbasi et al. 2008) and Auger (Abraham et al. 2008d) Collaborations, together with the Auger bounds on the UHE neutrino flux (Abraham et al. 2008b), on the photon fraction (Abraham et al. 2008a), and on the anisotropy toward the Galactic center (Aglietta et al. 2007; Aloisio & Tortorici 2008) are all consistent with this scenario. The next step is clearly the identification of the sources of UHECRs, an arena where anisotropy studies play a crucial role. Yet, the limited angular resolution of extensive air shower detectors and especially the deflections that charged particles suffer in astrophysical magnetic fields make the task highly nontrivial. This is especially troublesome given that the UHECRs chemical composition is unknown, that we lack a detailed knowledge of the Galactic magnetic field structure, and, above all, of the very magnitude and structure of extragalactic magnetic fields (EGMF) outside of cluster cores. These limitations—together with the small statistics available at present—suggest that, at least in an initial phase, charged particle astronomy may be limited to the inference on the statistical properties of UHECR sources rather than a detailed study of single accelerators.

In Cuoco et al. (2008b), we found that a global comparison of the two-point autocorrelation function of the data with one of the catalogs of potential sources is a powerful diagnostic tool: this observable is less sensitive to unknown deflections in magnetic fields than the cross-correlation function, while keeping a strong discriminating power among source candidates. In particular, the autocorrelation function of (sub-) classes of galaxies has different biases with respect to the large-scale structure (LSS) of matter. As a result, the best-fit value for the density n_s of different source classes may differ, especially if only one or a small range of angular scales is considered. Although the bias of different source classes differs maximally at small angular scales, we showed that the statistically most significant differences are at intermediate angular scales, where both the larger number of cosmic ray (CR) pairs and of galaxy pairs lead to relatively smaller error bars. Moreover, the autocorrelation function on larger angular scales becomes less dependent on possible deflections in the Galactic and extragalactic fields.

In this paper, we derive the number density of UHECR sources using the recently published arrival directions and energies of the 27 Auger events (Abraham et al. 2008c) with estimated energy E ⩾ 57 EeV, thereby complementing the study (Cuoco et al. 2008b) with a concrete example for a comparison of the global cumulative autocorrelation function of sources and UHECRs. Note that we showed in Cuoco et al. (2008b) that, even in an idealized case where systematics play no major role, roughly three times the number of "useful" events that can be extracted from Abraham et al. (2008c) are required to start distinguishing between different subclasses of sources. Thus, a study of the kind envisaged in Cuoco et al. (2008b) is unrealistic at present. Here, we restrict ourselves to more modest goals: (1) to compare the data against predictions of two toy model cases of uniformly distributed sources and of "normal galaxies" (i.e., sources that have the same distribution as the PSCz catalog (Saunders et al. 2000)) which we shall refer to with the two values for the label κ = {uni, LSS}, respectively. (2) To study the effect on the allowed range of n_s of a systematic error on the energy scale of the UHECR experiment. Note that preliminary results of the clustering of the Auger events has been presented in Mollerach (2007), but astrophysical implications have not been discussed there.

We review the statistical method that we use in Section 2, and apply it in Section 3 to the Auger data, providing some interpretation of the results. In Section 4, we discuss some limitations and caveats of the analysis. Finally, we summarize in Section 5.

The use of correlation functions is well suited to the study of over- and underdensities of nonuniformly distributed sources and of the resulting anisotropies in the radiation received from them. Since the number of CR events published by Auger is still small, we use in our analysis, following Kachelrieß & Semikoz (2006), the cumulative two-point autocorrelation function ${\mathcal C}(\vartheta)$ defined as

$\begin{equation} {\mathcal C}(\vartheta) = \sum _{i=2}^N\sum _{j=1}^{i-1}\Theta (\vartheta -\vartheta _{ij}) \,, \end{equation} \tag{ 1 }$

i.e., the number of pairs within the angular distance ϑ. Here, N is the number of CRs considered, ϑ_ij is the angular distance between events i and j, and Θ is the step function (with Θ(0) = 1).

This function is straightforward to compute for the data, and denoted then by ${\mathcal C}_\ast (\vartheta)$. For a specific model hypothesis X, a set of functions ${\mathcal C}_i(\vartheta |X)$ is obtained in the following way: sources with equal luminosities are distributed within a sphere of 180 h⁻¹ Mpc either uniformly or following the three-dimensional LSS as given by a smoothed version (Cuoco et al. 2006, 2008a) of the PSCz catalog (Saunders et al. 2000). For the LSS case (but not for the uniform case) sources and CR events within the PSCz mask are excluded, leaving 22 CR events. Note that the mask mostly overlaps with the Galactic plane and the bulge region, where larger deflections due to the Galactic magnetic field are expected: the mask is thus not only a catalog limitation, but also implements to some extent a physically motivated angular cut. Finally, each source k, at redshift z_k, is weighted by the factor

$\begin{equation} \frac{1}{z_k^2} \int _{E_i(E_{\rm cut},z_k)}^{\infty } E^{-s}\, {d}E \end{equation} \tag{ 2 }$

accounting for its redshift-dependent flux suppression and the CR energy losses. These are parameterized as a continuous process in term of the function E_i(E_f, z_k) which gives the initial injection energy E_i for a particle leaving the source at z_k and arriving at the Earth with energy E_f < E_i. Further details are given in Cuoco et al. (2006). The injection spectral index s is assumed to be the same for all the sources and equal to 2.0. The dependence on s is however weak as shown in more detail in Cuoco et al. (2006). This procedure defines the model, while a single random realization is obtained by choosing the observed number of events from the sources according to the source weights and the declination-dependent Auger experimental exposure.

The model thus directly depends only on n_s, and the choice between sources distributed uniformly or according to the PSCz catalog (of course, implicitly it also depends on the assumptions of sufficiently small magnetic field deflections). Additionally, the model depends via the weights of Equation (2) on the type of the primary particle, the energy spectrum of the sources, and the energy scale and cut E_cut. The latter dependence arises, because the energy scale of CR experiments has a relatively large systematic error that is difficult to determine. In particular, it has been argued (Teshima 2007) that the energy scale of Auger should be shifted up by 30%–40% to obtain agreement with the spectral shape predicted by e⁺e⁻ pair production (Berezinsky et al. 2006) and the CR flux measured by other experiments. Therefore, we use two different values for the energy cut, E_cut = 60 EeV, assuming that the energy calibration of Auger is correct and E_cut = 80 EeV inspired by the dip interpretation. We do not include in this work the finite energy resolution of the Auger experiment that is of order 20% in ΔE/E. A finite energy resolution would result in an effective decrease of the nominal energy cut due to the steeply falling CR spectrum and to a larger GZK horizon (Kachelrieß et al. 2007). Similarly, we do not account for the stochasticity of the energy loss in the photopion production process. Both of the effects are subdominant at the moment with the present low statistics while a more careful analysis will be needed in the future when more data are available. For a rough estimate of the influence of both effects on n_s one may compare how the 30% upward shift of E_cut from 60 to 80 EeV changes n_s. Finally, throughout this work we consider proton primaries, but note that the combination of nuclei with large deflections and few sources has been advocated too (Armengaud et al. 2005; Fargion 2008; Gorbunov et al. 2008). A further brief discussion of this point is postponed to Section 4.2.

The cumulative autocorrelation of the data ${\mathcal C}_\ast (\vartheta)$, which is a single, one-dimensional function, has now to be compared with the hypothesis X, for which we have various Monte Carlo realizations ${\mathcal C}_k(\vartheta |X)$, k = 1, ..., M (we use typically M ∼ 10⁵). A standard method to compare data and model is to use angular bins ϑ_i so that to substitute the continuous function ${\mathcal C}(\vartheta)$ with the discrete set of values ${\mathcal C}(\vartheta _i)$. The Monte Carlo realizations can then be used to calculate the marginalized probability distribution of each single ${\mathcal C}(\vartheta _i|X)$ or, if required, the joint probability distribution of two ${\mathcal C}$ variables ${\mathcal C}(\vartheta _i|X),\; {\mathcal C}(\vartheta _j|X)$, three ${\mathcal C}$ variables or more. In practice, to derive the best-fit value n_s as well as the goodness-of-fit for the chosen hypothesis X a possible way is to calculate the mean $\langle {\mathcal C}(\vartheta _i|X)\rangle$ and the variance σ_i per bin as well as the correlation matrix σ_ij and then to perform a χ² test. However, the difficulty to deal with such a high-dimensional probability space and the generally strong non-Gaussianity of the probability distribution make the χ² method clearly not optimal for the problem at hand.

The usual way to circumvent the problem is to use the Monte Carlo to calculate the chance probability to observe stronger clustering than in the data. Given the problem at hand we slightly generalize the method defining two functions: the chance probability to observe stronger (P₊(ϑ|X)) or weaker (P₋(ϑ|X)) clustering at the angular scale ϑ than in the data given by

$\begin{equation} P_{\pm }(\vartheta |X) = \frac{1}{M}\sum _{i=1}^M \Theta [\pm (C_i(\vartheta |X)-C_{\ast }(\vartheta))] \,. \end{equation} \tag{ 3 }$

The minimum of the chance probability can then be used as a global estimator for the agreement of the hypothesis with the data. In our case with a scan of P_±(ϑ|X) over ϑ, we obtain the two minima P₋(X) for the minimum of P₋(ϑ|X) and P₊(X) for the minimum of P₊(ϑ|X), respectively. The simplest way to combine these two estimators is then to use the product P(X) = P₋(X)P₊(X).

However, a drawback of this method is that neither the quantity P(X) nor the single P_±(X) are truly probabilities. The scan over the angular scales ϑ, in fact, introduces a bias which need to be corrected with the use of a penalty factor (Tinyakov & Tkachev 2001, 2004; Finley & Westerhoff 2004). More precisely, to obtain correct probabilities the identical procedure as described above needs to be performed with many simulated data sets. Counting how often smaller values of P_±(X) and P(X) are obtained by chance with respect to the case of the data provides true penalized chance probabilities p₊(X), p₋(X), and p(X).

The use of the chance probability tool has often generated some confusion in the past. An emblematic case is the significance of the small-scale clustering in the AGASA data for which very different estimates ranging from ∼10⁻⁶ to ∼10⁻² have been reported in various studies (see for example Takeda et al. 1999) depending on the use or not of the penalty and on different a priori choices of the angular and energy scale of reference (see Finley & Westerhoff 2004 for a detailed account). The effect of the scan can then be quite relevant and a major point in the following is that the effect of the penalty is correctly taken into account when quoting the constraints on n_s and the constraints are further compared with the case of an a priori choice of the relevant angular scale.

Note, anyway, that the penalty calculation can be avoided if a particular angular scale is chosen a priori and the values of the chance probability at this scale are employed. However, the scan over all angles avoids possible bias, in contrast to the choice of a single angular scale, which introduces some theoretical prejudice even if this choice may be physically motivated. In the case at hand, it is unclear if this should be dominated by the ∼1° angular resolution of the detector, as in Takami & Sato (2009) which implicitly assumes negligible magnetic field deflections; by a ∼3° scale, as suggested by the cross-correlation with active galactic nuclei (AGNs) revealed by Auger (Abraham et al. 2007);⁵ or yet some other scale, as the 6° separation considered in Abraham et al. (2008c). To emphasize this dependence, we summarize the results of this kind of "single-bin" analyses in Table 1 and compare them with the ones obtained with the global method, as reported in the last column (see the next section). Note, in particular, how the 6° bin chosen in Abraham et al. (2008c) systematically leads to an overly stringent bound.

Table 1. The Estimated Number Density of Sources (at 95% Confidence Level) Under Different Assumptions, Using Only the First Bin Information with Different Sizes, and Compared with the Global Method

ns(10−4 Mpc−3)	ϑ1
	6°	$3\sqrt{2}^\circ$	3°	1°	Global
LSS (80)	1.3+5.7−1.0	2.0+8.0−1.6	2.0+∞−1.4	5.0+∞−4.2	1.3+100−0.8
Uniform (80)	0.8+0.8−0.5	1.2+0.8−0.8	2.5+∞−1.8	5.0+∞−4.2	1.4+1.4−0.7
LSS (60)	0.3+1.7−0.27	0.3+2.7−0.18	0.7+100−0.5	1.5+100−1.25	0.8+19−0.6
Uniform (60)	0.2+1.0−0.12	0.3+1.7−0.2	0.8+70−0.63	1.0+∞−0.8	0.5+0.5−0.2

Download table as:  ASCIITypeset image

In Figure 1, we report the results for the quantities p_i(X) defined above, where X ≡ {n_s, E_cut, κ}, in our case the latter two being two-valued discrete variables. Because the number of CRs usable for this analysis is still small, all four hypotheses are compatible with the data at the 2σ level for some range of n_s values. Yet, several interesting conclusions can already be drawn. The best fit is achieved for sources following the PSCz distribution and a source density n_s ≃ 1 × 10⁻⁴ Mpc⁻³. Also, independently of E_cut, both the penalized probability and the range of n_s with compatibility at 95% CL are larger for the LSS model than for the uniform case. Reversing the argument, as can be read more quantitatively in Table 1, we can see that the constraints are generally stronger for the uniform cases with respect to the LSS ones, but this is achieved only at the expense of a worse general best fit. The fact that the LSS models fit better than the data is not surprising: most of the Auger events are aligned along the local overdensity known as the supergalactic plane which is suitably reproduced with the use of the PSCz catalog within our LSS scenario.

Zoom In Zoom Out Reset image size

The case of a uniform distribution of "infinitely many" sources (n_s → ∞ each with an infinitesimal luminosity) is excluded for both energy cuts at the 95% CL: the upper bound is n_s ≲ (1–3) × 10⁻⁴ Mpc⁻³. This is another way to say that the Auger data are inconsistent with a structureless UHECR sky, independently of the use of a catalog and of a pre-determined angular scale for the search. This is, in our opinion, an important milestone in the development of UHECR astronomy. While the best-fit point for n_s is approximately a factor 10 higher than that found in earlier studies using AGASA data above E_cut > 40 EeV (in the AGASA energy scale; Yoshiguchi et al. 2003; Blasi & De Marco 2004; Kachelrieß & Semikoz 2005), the shape of the chance probability p(n_s) agrees: for low values of n_s, p(n_s) is a steeply decreasing function of n_s, since the probability to observe multiplets from the same source increases fast. In particular, the radius within which 70% of all observed UHECRs with energy above E_cut = 80 EeV are produced is R ≈ 60 Mpc (Cuoco et al. 2006).⁶ As a result, the number of sources within this radius becomes less than the number of observed CR events for densities smaller than n_s ≈ 10⁻⁵. Such a scenario would require large deflections (and probably nuclei primary) and thus contradicts our assumptions. On the other hand, p(n_s) decreases relatively slowly for high densities and only weak constraints can be obtained with the current data set for the maximally allowed value of n_s. Since both an increase of E_cut and of the bias of the sources lead to a decrease of the effective number of sources inside the GZK volume, large values of n_s have the strongest constraint in the case of uniformly distributed sources and E_cut = 60 EeV (left, top panel) and weakest for sources following the LSS and E_cut = 80 EeV (right, bottom panel).

In Figure 2, we show the model autocorrelation function with 1σ and 3σ shaded regions for the best-fit uniform and LSS model for E_cut = 80 EeV, both corresponding to n_s ≈ 1.4 × 10⁻⁴ Mpc⁻³, together with the data. At small angular scales, ϑ ≲ 3°, the data show a deficit of clusters compared to the expectation for the best-fit density from the global analysis. This "tension" is qualitatively present in most models fitting the data. Within our assumptions, this deficit is explained in a natural way by (relative) deflections of this size in magnetic fields. This value is comparable with the 3.2° found in Abraham et al. (2008c) that optimizes the correlations of the same data set with AGNs. The absence of small-scale clusters is also responsible for the shift in the best-fit value for n_s compared to old analyses using the AGASA data. The result thus shows that, intriguingly, also the autocorrelation function can be employed as a sensitive tool for magnetic field studies. Clearly, however, a more detailed study needs to be complemented with a model of the intervening magnetic fields while, likely, more statistics is required to derive significant constraints.

Zoom In Zoom Out Reset image size

Our result is in contrast with that of Takami & Sato (2009) which instead reports a small-scale clustering in the Auger data. Notice however that, with respect to our work, the authors of Takami & Sato (2009), although including an explicit treatment of galactic and extragalactic magnetic fields, make use of the noncumulative autocorrelation functions and do not take into account penalties for the scan over the angular scale. Their claim of small-scale clustering within 1° is thus likely to disappear if the angular scan is taken into account and the comparison is made properly with respect to a model effectively fitting the data. Thus, as already noted in Harari (2004), the interpretation of small-scale clustering could be misleading if it is not defined with respect to a proper model of the distribution of sources. Apart from this point, however, the constraints we obtain on n_s are roughly in agreement with Takami & Sato (2009). This is mainly due to the very low statistics available at present which is still not sensitive to the exact analysis method employed. We stress however, as noted in Cuoco et al. (2008b), that already with a statistic of ∼3 times the present one a formally correct analysis will become crucial to avoid biased results.

Figure 2 also clarifies why the LSS case gives a better fit to the data: as can be seen, the LSS best-fit model fits nicely to the data, basically within 1σ over all the angular scales, while the best fit case for uniform sources shows a ∼2σ deficit in the broad range 4°–30°. A lower n_s (i.e., a higher clustering) would not help because it would give much more pairs than the data in the 1°–4° range. Thus, at the end, even with the best possible compromise the agreement with the data is only at the 20% level for the best fit (or, equivalently, the uniform model is excluded at the 80% CL).

We summarize in Table 1 the list of best-fit n_s with 95% error bars for the four cases considered and for different choices of the angular scale or for the global autocorrelation analysis. The crucial point to note here is that the derived n_s intervals are sensibly biased with respect to each other for different choices of the a priori angular scale. Thus, the choice of the angular scale crucially affects the result and should be avoided unless the given scale has some strong physical motivation. This problem is of course avoided employing a global comparison. The choice of a single angular scale also affects the error bars which can be both larger or smaller with respect to the global case. This can be easily understood from Figure 3 where we can see that the choice 4°–15° is optimal because excess of clustering is observed for low n_s while a deficit is present for high n_s giving thus the tightest constraints. Again, however, a choice in this range is not a priori motivated so that from the global analysis we get somewhat larger error bars properly taking into account all the angular scales.

Zoom In Zoom Out Reset image size

Although it is difficult to draw strong conclusions on the nature of the UHECR sources at this moment, it is worth noting that X-ray-selected AGNs with X-ray luminosity L > 10⁴³ erg s⁻¹ naturally fall in the range n_AGN = (1–5) × 10⁻⁵ Mpc⁻³ (Steffen et al. 2003). At the same time, for AGNs with densities in the range (10⁻⁵–10⁻⁴) Mpc⁻³, the required luminosity is of the order ${\cal L}/n_{\rm AGN}\sim (10^{40}\hbox{--} 10^{41})$ erg s⁻¹ in UHECRs above ∼ 60 EeV. So, simple energetics arguments are consistent with the inferred values for the density; on the other hand, acceleration efficiency (see, e.g., Ptitsyna & Troitsky 2008) tends to favor some subclasses of objects and thus a lower inferred density of sources, yet typically within the 95% CL range inferred for n_s. Another effect which might reconcile the small tension is that the sources can be bursting/transient or beamed. In this case, the luminosity requirement is softened (see, e.g., the model proposed by Farrar & Gruzinov 2009). However, the effectively visible number density of sources decreases too, unless some isotropization takes place after acceleration but close to the source (see Sigl 2008, for more details). Of course, these are very preliminary considerations based on a limited number of data. For example, these arguments might be significantly modified when accounting for a more realistic luminosity function. The effects of extragalactic magnetic fields and the consequent UHECRs time delays can also be quite relevant, although at present our knowledge of EGMFs is affected by large uncertainties and needs to be better understood (see Sigl et al. 2003, 2004; Dolag et al. 2004; Kotera & Lemoine 2008; Ryu et al. 2008). Also, in the presence of a heavy-nuclei component and/or of extragalactic magnetic fields, a significant fraction of events might be associated with the nearest AGN Cen A (see, e.g., Gorbunov et al. 2008). Interestingly, signatures in the gamma and neutrino bands are expected to help disentangling among many scenarios (see Sigl 2008; Cuoco & Hannestad 2008; Kachelrieß et al. 2009; Becker & Biermann 2009; Halzen & O'Murchadha 2008; Hardcastle et al. 2008), enlarging the realm of multi-messenger astronomy.

4.1. The Role of the Energy Cut

The energy cut used by the PAO to produce the sample used in the present analysis was chosen in order to maximize the cross-correlation signal, see Abraham et al. (2008b, 2008c) for details. One may wonder how this selection affects the conclusions of our present work. In principle, the optimal energy cut for the autocorrelation signal and the cross-correlation signal with a given catalog differ from one another, although in the case of a common physical origin one does expect some correlation between them. In particular, it seems reasonable that the optimal cut for a global autocorrelation might reside at a lower energy, since the small-scale displacement from putative sources is not as relevant to the signal as for cross-correlations, and the larger statistics helps. Lacking a direct access to the data, it is hard to estimate quantitatively how large a bias is introduced by focusing on the sample of public available events. From Figure 2 presented in Mollerach (2007), however, one can draw two qualitative conclusions: (1) that a correlation between the two optimal cuts in energy is indeed present; (2) that a slightly different cut, in the range 40 ≲ E_cut/EeV ≲ 60, should still lead to an appreciable clustering of the events.

We checked that this is indeed the case by performing the same analysis as before, but now adding a further scan over the range of accessible energy cuts, i.e., any value E > 57 EeV. We plot the results in Figure 4.⁷ One can note that in all cases the best fit improves and the constraints on n_s worsen, as expected given the further penalty due to the energy scan. Yet, most qualitative features described in the previous section stay the same: for example, the LSS model is still preferred over a uniform one. It should be also said that if the true minimum of the chance probability is below E = 57 EeV and thus not included in the scan, then these constraints are "overpenalized" and thus looser than necessary. At the moment, it is impossible to draw more quantitative conclusions, since our scan suggests that it is likely that the optimal cut for the autocorrelation function is below E = 57 EeV, a range for which the events are not publicly available.

Zoom In Zoom Out Reset image size

4.2. Assumptions on the Chemical Composition

We used the first UHE data released by Auger to perform a global autocorrelation study of UHECR arrival directions, assuming proton primaries. The major advantage of our tool is that no biases are introduced by the a priori choice of a single angular scale. The main observable we have focused on is the number density n_s of ultrahigh energy cosmic ray (UHECR) sources. While the global analysis does not bias n_s by what is the theoretical prior of the "relevant" angular scale, still it is important to establish how the extraction of the allowed range of n_s from the data depends on a number of other effects, an issue often overlooked in the literature. In particular, here we discussed the systematic energy scale uncertainty and the bias of UHECR sources with respect to LSSs. As a first attempt to extract some information from the data, we compared four hypotheses: a structured universe (following the PSCz catalog) and an isotropic case, each for two possible values for the absolute energy scale of the PAO experiment. The density n_s is the free continuous parameter in terms of which constraints have been analyzed.

Not surprisingly, we find that the number of CRs usable for this analysis is still small and all four hypotheses are compatible with the data at the 2σ level for some range of n_s. Yet, several interesting observations can be tentatively drawn. The best fit is achieved for sources following the matter tracer distribution, and a source density n_s ≃ 1 × 10⁻⁴ Mpc⁻³. It is interesting to note that the data show some preference (although not significant, yet) for a structured universe: other recent statistical studies, like Kashti & Waxman (2008) and Koers & Tinyakov (2009), qualitatively agree in that respect. Also, there is an indication that the case of a uniform distribution of "infinitely many" sources (n_s → ∞ each with an infinitesimal luminosity) is excluded for both energy cuts at the 95% CL: the upper bound is n_s ≲ (1–3) × 10⁻⁴ Mpc⁻³. This is another way to say that early Auger data suggest that data are poorly consistent with a structureless UHECR sky, independently of the use of a catalog and of a pre-determined angular scale for the search.

Compared to a benchmark number density of proton sources n_s ≃ 1 × 10⁻⁴ Mpc⁻³, a factor ∼2 lower densities are preferred if the current Auger energy scale is correct, a factor ∼2 higher value if it is underestimated as required by the dip model. Including the finite energy resolution of the Auger experiment into the analysis will further reduce the best-fit value for n_s. The width of the allowed region is dominated at present by the statistical error due to the small number of events. Nominally, approximately a three times larger sample is needed to reduce the Poisson error below the typical differences between source candidates. Once that level of statistics is reached, other effects will provide the main source of error, a major one being the systematics on the energy scale as we illustrated here. In the future, other effects such as the systematic and statistical errors in the energy determination of UHECR events and the stochastic nature of the photopion process need to be included to correctly determine the best-fit value of n_s. Even considering the above limitations, preliminary conclusions are the following. First, the fit generally improves for sources following the LSS compared to uniformly distributed sources; qualitatively, for AGNs which are known to be even more overdensity biased than the LSS (see, e.g., discussion in Cuoco et al. 2008b), the agreement is expected to improve even further. In particular, the case of an uniform distribution of "infinitely many" sources, i.e., a structureless UHECR sky, is excluded at 99% CL. Second, the absence of clustering on scales smaller than a few degrees is most easily understood as the effect of magnetic smearing of comparable size. This, intriguingly, suggests that autocorrelation studies can be employed as a complementary tool to study galactic (and extragalactic) magnetic fields. Conclusions about the source scenario are complicated by the limited knowledge of the EGMFs and the presence of beaming and/or bursting effects which are difficult to disentangle with the use of the autocorrelation alone. Hence, the true source density could be somewhat lower than the best-fit value found: the energy scale error alone shifts the best-fit value by a factor 2 or 3; luminosity function effects are likely relevant too. The scarce statistics at the moment is a serious limiting factor to constrain more realistic models, but with a greater exposure of currently existing instruments and the multi-messenger combination of gamma-ray and neutrino data, some of these issues will probably be addressed and solved in the near future.

For a significant contamination of nuclei in the sample, a much more complicated analysis is needed, since many more variables enter the game. Qualitatively, one can expect that although future data might allow to disentangle a proton-dominated sample from a more complicated picture, inferring source properties and disentangling them from magnetic effects (in a few words, performing UHECR astronomy) should wait for a major jump in exposure, perhaps beyond the capabilities of currently planned instruments. Similar considerations apply unfortunately to the constraints to the cross sections of UHECRs as well.

We are grateful to Günter Sigl for useful discussions. M.K. thanks the Max-Planck-Institut für Physik in Munich for hospitality and support. P.S. is supported by the US Department of Energy and by NASA grant NAG5-10842. Fermilab is operated by Fermi Research Alliance, LLC under contract no. DE-AC02-07CH11359 with the United States Department of Energy.