Testing large-scale (an)isotropy of ultra-high energy cosmic rays
We present a simple yet powerful method to test models of cosmic-ray (CR) origin using the distribution of CR arrival directions. The method is statistically unambiguous in the sense that it is binless and does not invoke scanning over unknown parameters, and general in the sense that it can be applied to any model that predicts a continuous distribution of CRs over the sky. We show that it provides a powerful discrimination between an isotropic distribution and predictions from the “matter tracer” model, a benchmark model that assumes small CR deflections and a continuous distribution of sources tracing the distribution of matter in the Universe. Our method is competitive or superior in statistical power to existing methods, and is especially sensitive in the case of relatively few high-energy events. Applying the method to the present data we find that neither an isotropic distribution nor the matter tracer model can be excluded. Based on estimates of its statistical power, we expect that the proposed test will lead to meaningful constraints on models of CR origin with the data that will be accumulated within the next few years by the Pierre Auger Observatory and the Telescope Array.
The recent observation [1, 2] of a strong suppression in the energy spectrum of ultra-high energy cosmic rays (UHECRs) provides compelling evidence that the bulk of UHECRs come from relatively close (within Mpc) extragalactic sources. It seems inevitable that the distribution of cosmic-ray (CR) sources traces the distribution of matter on these length scales. One should therefore expect anisotropy in the UHECR flux due to the non-uniform distribution of matter in the nearby Universe. Although other experiments have not detected deviations from isotropy at the highest energies, the Pierre Auger Observatory (PAO) has recently reported anisotropy in the flux of UHECRs above 57 EeV at a confidence level of more than 99% [3, 4].
Exactly how CR anisotropy manifests itself depends crucially on two factors: the deflections of CRs in magnetic fields and the density of sources. The deflections of UHECRs in Galactic and extragalactic magnetic fields may vary from a few to a few tens of degrees, the uncertainty being due to both unknown parameters of the magnetic fields and the uncertain chemical composition of CRs at the highest energies. The extragalactic magnetic fields are usually assumed to be small ( G) in voids, and larger in clusters of galaxies and filaments. Simulations show that deflections in the extragalactic magnetic fields may be small for protons for most directions on the sky  (see, however, Ref. ). For iron they may reach several tens of degrees. The deflections in the Galactic magnetic field may also vary from a few to a few tens of degrees, depending on the chemical composition of UHECRs. However, even the largest deflections are unlikely to completely wash out the anisotropy created by the non-uniform source distribution.
Although the CR source density affects the anisotropy signature, deviations from isotropy are expected for any density: If the number of sources of observable CRs is large, anisotropies in the CR flux arise from the fact that the distribution of sources is inhomogeneous because it traces that of visible matter. If, on the other hand, there are few sources, they will produce anisotropy just because of that. Anisotropy is therefore expected for any source density.
Different assumptions on deflections and source densities lead to different predictions for the UHECR anisotropy. The case of small deflections is of particular interest because in this case the uncertainties related to magnetic fields are not present. Under the assumption of small deflections, the existing data imply that the number of observable sources above EeV is large, namely several hundred or more [4, 7, 8]. This situation may be approximated by a continuous distribution of sources following the matter distribution in the Universe. This model, which we refer to as the “matter tracer” model, predicts that local matter overdensities such as nearby galaxy superclusters should be visible in UHECRs. The matter tracer model may serve as a benchmark in anisotropy studies, as its main ingredients – small deflections and a large number of observable sources – represent a limiting case of many realistic situations. Apart from details on particle acceleration in CR sources, the model has no free parameters: for a given injection spectrum the distribution of UHECR over the sky is predicted uniquely. By comparing this prediction to observations, one may test to which extent the underlying assumptions — small deflections in the magnetic fields, in the first place — are correct.
The prediction of a continuous distribution of CRs over the sky calls for special statistical tests. In this paper we propose a new statistical test suitable in this situation. Several other methods have been discussed in the literature [6, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Our method is competitive or superior in statistical power to the existing methods (see Sect. 4 for a comparison in power between different methods). It is especially sensitive in the case of a few high-energy events, which is a regime of physical interest.
The idea of the method is illustrated in Fig. 1, which shows the actual flux distribution above 40 EeV predicted in the matter tracer model together with 15 fictitious, uniformly distributed CR arrival directions. (We discuss how this figure is obtained in Sect. 2.4). The contours of equal flux are drawn in such a way that they divide the sky into 5 areas, each containing an equal share of the total flux (1/5 in this case). This implies that events produced according to the matter tracer model should be distributed roughly equally over the five bands. On the contrary, if the sources are distributed in a different way, one generally expects deviations from equipartition. For instance, in the case of an isotropic CR flux most events will fall in the region with the smallest flux (the lightest band in Fig. 1) just because this region covers the largest area. Comparing the distribution of observed events over the five bands to a model prediction then allows one to test the underlying model. This is illustrated in the inset of Fig. 1 with a histogram that shows the number of events in every band. For this specific example one would reject the matter tracer model as the distribution of events over the five bands deviates significantly from equipartition.
The method described above can be formulated in a binless way. This is an important advantage because binning introduces ambiguity which complicates the statistical analysis. Our numerical method starts by generating a (continuous) skymap of expected flux. We then extract flux values at the positions of data events, thus obtaining a set of flux values which we refer to as the “data set”. Then a large number of simulated events is generated by Monte Carlo event generation from the model that is to be tested (e.g., an isotropic distribution or the matter tracer model). The flux values at positions of simulated events form the (Monte Carlo) “reference set”. Then a statistical test is performed to determine whether the data and reference sets are drawn from the same distribution. This final step can be done with the well-known Kolmogorov-Smirnov (KS) test. In this work we also propose a modified version of the KS test, which we show to be more sensitive than the standard KS test for the case at hand. Both the KS test and its modification do not require binning, rendering the whole procedure binless. Furthermore, the method does not invoke scanning over unknown parameters and hence avoids the use of penalty factors to estimate statistical significances.
Although the main purpose of this paper is to introduce our statistical test and analyze its potential to infer properties of UHECR sources from future data, we have also applied it to the publicly available data from the Akeno Giant Air Shower Array (AGASA) and PAO to test the matter tracer model and an isotropic CR distribution. As neither of the two hypotheses can be ruled out, the present data give inconclusive results which do not significantly add to the ongoing debate on the sources of UHECRs. We find, however, that the proposed test offers a powerful discrimination between source models with more CR events, and we believe it is timely to present the method before such data become available.
The rest of this paper is organized as follows. In Sect. 2 we discuss our modeling of CR arrival directions in the matter tracer model. This includes subsections on the modeling of the source distribution, CR propagation, and on the Monte Carlo event generation procedure. Sect. 3 regards anisotropy tests. In this section we discuss the test statistics considered in this work, and recapitulate some essential concepts in statistical hypothesis testing. In Sect. 4 we compare the statistical power of the anisotropy test proposed in this work to other tests and we consider its dependence on model parameters. We apply our test to data from AGASA and PAO in Sect. 5. We discuss our work in Sect. 6 and conclude in Sect. 7.
2 Modeling of CR arrival directions
In this work we will focus on two different models for the distribution of CRs. The first is the matter tracer model, which assumes that CR sources trace the distribution of matter and are numerous. The second model (the “isotropy” model) predicts an isotropic distribution of CR events over the sky, barring possible non-trivial detector exposure. In this section we discuss our modeling of the CR arrival directions in the matter tracer model (the modeling of arrival directions in the isotropy model being trivial).
2.1 Source distribution; the KKKST catalog
To estimate the distribution of CR arrival directions in the matter tracer model we must first model the distribution of their sources. The matter tracer model presumes that CR sources trace the large-scale distribution of matter in the Universe. To our knowledge the most accurate way to model the large-scale structure is by using galaxy catalogs obtained from extensive surveys. There is however some ambiguity in this procedure because different catalogs lead to different reconstructed matter density distributions. The reasons for this are twofold. First, different types of galaxies exhibit somewhat different clustering properties so that selection effects will affect the reconstructed matter distribution. Second, each of the available catalogs has certain drawbacks that limit the reconstruction accuracy. These include incompleteness of the catalog, limited sky coverage, and uncertainties in redshift determination.
For the present study we have considered two different galaxy catalogs: the PSCz catalog  and a catalog that was recently compiled by Kalashev et al. , hereafter referred to as the KKKST catalog. The PSCz catalog contains galaxies up to 4 Gpc. It has the advantage of accurate redshift determinations, but it is is incomplete and it does not cover all regions of the sky. Catalog incompleteness may be taken into account by a selection function; however this necessarily invokes some extrapolation. The incomplete sky coverage is also a disadvantage because it necessitates either interpolating the matter density over the left-out regions or a specialized treatment to include only a specific part of the sky.
The KKKST catalog, compiled from the 2MASS XSC catalog  and the (HYPER)LEDA database , contains galaxies within 270 Mpc. It is a complete sample up to 270 Mpc for galactic latitude . The 2MASS XSC catalog is at present the most complete all-sky galaxy catalog. Radial velocity measurements are however unavailable for the bulk of the galaxies, so that one has to rely on photometric redshift estimates (assuming galaxies are standard candles in the infrared).111The XSCz catalog (see http://web.ipac.caltech.edu/staff/jarrett/XSCz/) contains radial velocity measurements for 27% of the 2MASS galaxies. It is however not (yet) publicly available. As discussed in Ref. , photometric redshift estimates lead to large uncertainties in individual distance estimates but yield correct estimates when averaged over a sufficiently large number of galaxies. For galaxy groups and clusters, the accuracy is estimated to be within % . Although the 2MASS XSC galaxy catalog is expected to lead to accurate reconstructed matter densities at large distances, where there are many galaxies, it may not be reliable at smaller distances. For this reason the KKKST catalog uses the (HYPER)LEDA database, which contains measured radial velocities, at distances below 30 Mpc.
In this work we will mostly use the KKKST catalog to reconstruct the local matter density in the Universe, choosing inaccuracy in redshift determination over inaccuracies associated with catalog incompleteness. (In Sect. 4 we will also recompute some of our results using the PSCz catalog to study the effect of the assumed galaxy catalog on the power of the statistical tests proposed in this work.) We find the inaccuracies due to photometric redshift estimates acceptable, given the other uncertainties we are faced with. In particular, the uncertainties in CR energies due to experimental resolution are of the order 20%, which translates to an uncertainty in viewing distance (and hence on the distance scale) %. To assess the differences between the PSCz and the KKKST catalogs, we have compared the CR fluxes predicted by the matter tracer model with the matter density in the Universe reconstructed from both catalogs. We find, reassuringly, that the flux distributions over the sky are similar in the sense that the same physical structures may be recognized (see also below), although some are brighter in one catalog compared to the other. The most important difference between the two is that the KKKST catalog gives rise to stronger contrasts between under- and overdense regions. This also implies a stronger contrast in expected CR fluxes and hence in larger differences between predictions from the matter tracer model and an isotropic CR distribution. Whether or not the true CR sources exhibit these relatively large contrasts is at present unknown. Our approach is to adopt the KKKST catalog, and hence the relatively large density contrasts, at face value and let current and future data decide whether or not the predictions derived from it are accurate.
The KKKST catalog contains five objects that are closer than 5 Mpc. We remove these from the catalog because they do not represent a proper statistical sample of the local structure, while they would outshine all other sources due to their proximity.222This notably excludes Centaurus A, with a distance of Mpc the nearest active galaxy, from our analysis. The PAO data raise the possibility that this galaxy is a source of UHECRs [3, 24, 25, 26, 27, 28, 29]. Other methods (e.g., small-scale anisotropy tests) are more appropriate to test this possibility.
The distribution of galaxies as a function of distance is shown in figure 2. As can be seen in the figure, the distribution is consistent with a uniform density at distance larger than a few ten Mpc, which reflects the completeness of the sample. At smaller distances there is an overdensity which we attribute to galaxies in the local supercluster. The global density of galaxies in the catalog is Mpc.
The absence of data at large distances and small galactic latitude (close to the galactic plane) necessitates extrapolation. At distances beyond 270 Mpc, we will assume that sources are distributed uniformly. We exclude the region within of the galactic plane from our analysis. To compensate for artificial boundary effects due to this cut, we extrapolate the catalog data in a straightforward manner when computing average source densities just outside the excluded region: for every galaxy with coordinates and , where , we add a mirror galaxy with coordinates and in the galactic strip. Here the plus (minus) sign corresponds to ().
In Fig. 3 we show the expected integral CR flux above 40 EeV in the matter tracer model, using the KKKST catalog to reconstruct the local matter density in the Universe, together with the positions of some prominent groups and (super)clusters of galaxies. (We discuss how this figure is obtained in Sect. 2.4). In this figure, as in all figures in this paper showing model CR fluxes on the sky, the gray bands are chosen such that they contain 1/5 of the model flux each, with darker bands indicating larger flux. (We would like to stress that this division in a discrete number of bands is for presentation only.) Both far-away superclusters and close-by galaxy groups can be recognized in the figure, the most prominent overdense region extending between the local (Virgo) and Centaurus superclusters.
2.2 Source density
The CR source density is unknown, but bounded from below by the scarcity of observed doublets and triplets, and also by the fact that the sources of the observed UHECRs cannot be very far away because of CR energy loss in the Universe. These considerations lead to a (rough) lower bound Mpc .
The CR source density affects the distribution of CR arrival directions through event clustering: if there are fewer sources, one expects more doublets, triplets, etc. In the limit of very many sources the CR distribution becomes independent of source density, as the probability to observe more than one CR from any source is small regardless of the exact value of the source density. We will refer to this situation as the many-source regime.
The matter tracer model considered in this work is built on the assumption of a large number of observable sources. Since we model the distribution of CR sources from the KKKST catalog, we should verify that the number density of our model sources is sufficiently large for the many-source regime to apply (if it were too small we would risk introducing spurious event clustering). Assuming equal source luminosities, the number of CR sources may be estimated from the number of singlets and doublets :
Demanding, somewhat arbitrarily, that we find that . This translates to a typical source density , above which the many-source regime applies:
where is the total number of CR events and denotes the maximum viewing distance (which is limited due to interactions between CRs and the CMB). For realistic numbers of events the source density in the KKKST catalog and the many-source regime applies. Hence we can model the distribution of CR arrival directions in the matter tracer model from the KKKST catalog.
In Sect. 4 we will be interested in the relation between source density and power of the statistical test proposed in this work. As a matter of convenience we parameterize the source density by the density fraction
The lower bound on the CR source density then implies that .
2.3 Cosmic-ray propagation
In this section we consider modifications to the CR spectrum due to redshift and interactions with the cosmic microwave background (CMB). Throughout this work we will assume that CRs are protons, the composition of CRs at the highest energies being unclear [30, 31, 32]. CR energy loss in the Universe may be described by the following differential equation (see, e.g., Ref. )
where is the proton energy, denotes the source redshift, is defined below, and denotes the inverse energy loss time (at present epoch) due to CR-CMB interactions. We include CR energy loss due to photopion production and electron-positron pair production. For these processes we use an approximation for , valid for energies above eV, which may be found in the appendix of Ref. . This simplified treatment does not account for the stochastic nature of CMB-CR interactions nor for the limited energy resolution in the detector (both effects lead to underestimation of the mean free path; see, e.g., Ref. ), and it neglects interactions with the extragalactic background light (which overestimates the mean free path; see, e.g., Ref. ). In the light of other limitations, notably the scarcity of data at the highest energies, we consider these simplifications justified in modeling the arrival directions of CRs. We note that these simplifications do not affect the statistical test in itself.
Within the CDM concordance model, the redshift and proper distance are related as follows:
where is the speed of light, denotes the present Hubble rate, and () is the present density of matter (vacuum energy) in units of the critical density. We use the concordance model values km s Mpc, , and .
Evolving eq. (4) forward in time (backward in redshift) yields the energy of a proton on its trajectory from source to Earth. In figure 4 we show the CR energy as a function of traversed distance for different injection (source) energies . The figure clearly shows the strong energy loss of ultra-high-energy protons that gives rise to the GZK cutoff in the CR energy spectrum.
The physical quantities relevant to the present analysis are easily expressed in terms of a function that gives the CR energy at the source as a function of observed energy and source distance . We determine this function by numerically evolving eq. (4) backward in time (forward in redshift). In terms of the integral flux from a single source at distance is:
where is the integral luminosity above of the source. Summing up the contributions of all sources, the diffuse CR flux may be expressed as follows:
where is the present source density, is the maximum redshift (we use unless otherwise stated), and parameterizes source evolution (no evolution corresponds to ). In this equation we assume, as we do throughout this work, identical sources.
2.4 Cosmic-ray flux distribution in the matter tracer model
In this section we express the CR flux from a direction on the sky in the matter tracer model using the KKKST catalog and the CR propagation model discussed in the previous subsections. We are interested in the integral flux above some threshold energy . As the catalog is limited to a catalog distance , we first separate the total flux into two parts:
where denotes the contribution that arises from sources within the catalog distance , denotes the contribution from sources at larger distance (which is isotropic by assumption), and and are galactic coordinates. We parameterize the relative strength of the two components as follows:
which defines as the fraction of the total CR flux above energy produced by sources within distance . We compute this fraction from eq. (7) by integrating up to redshift corresponding to and dividing the result by eq. (7) integrated to . This fraction is shown in the top panel of Fig. 5 for three values of . In producing this figure we have assumed an injection spectrum with extending to energies much larger than 100 EeV, and no source evolution. As may be seen from the figure, the fraction of sources beyond 270 Mpc contributing to the observed CR flux is essentially zero above EeV. In the bottom panel of Fig. 5 we show the effective CR viewing distance as a function of CR energy. The three lines show the distances within which 50%, 90% and 99.9% of the total flux are produced.
We model the contribution from sources within the KKKST catalog to the diffuse CR flux as follows:
where enumerates the galaxies in the catalog, is the integral flux above from galaxy (computed with eq. (6)), and
is a weighting factor that we insert to introduce angular smearing, being the angle between galaxy and the line of sight and being the angular smearing scale. The angular smearing procedure is introduced to avoid unphysical fluctuations due to the use of a catalog of point sources, and to account for CR deflections in the intergalactic magnetic fields and for limited angular resolution in the detector. We conservatively adopt in this work (see Ref.  for a detailed discussion).
In figure 6 we show the relative strength of the CR flux on the sky in the matter tracer model. Here we have used the KKKST catalog to reconstruct the local matter density in the Universe. The grayscale indicates the value of the quantity
where overlined quantities are averaged over the sky. In this figure we consider three different threshold energies: 40 EeV, 60 EeV, and 100 EeV. The gray bands are chosen such that each band contains 1/5 of the total flux, with a darker gray indicating a larger flux. A comparison between the different panels shows that the contrast in the expected flux strengths increases strongly with increasing CR energy. This is due to the fact that the nearby structure becomes more prominent with decreasing viewing distance.
2.5 Event generation
Here we discuss the Monte Carlo procedure employed in this work to generate CR arrival directions (“events”). The procedure relies on our modeling of the galaxy distribution (discussed in Sect. 2.1), the density of sources (Sect. 2.2), CR energy loss (Sect. 2.3), and on the detector exposure.
Event generation in the isotropy model is trivial: we determine arrival directions on the sky such that their distribution is proportional to the detector exposure.
In the matter tracer model CR sources are selected from the KKKST galaxy catalog. Here the situation is complicated by the fact that the catalog does not include sources beyond 270 Mpc. In case the maximum viewing distance is smaller than this (which is the case for CR energies above EeV), we proceed as follows. First we select a galaxy from the KKKST catalog, where the probability to select galaxy is
being the relative detector exposure in the direction of the event and the expected flux from galaxy above energy (cf. eq. (6)). Once a galaxy is selected, we generate an arrival direction around the direction of the galaxy such that the distribution of the angle between galaxy and arrival direction is
where is the smearing angle (see Sect. 2.4).
If the maximum viewing distance is beyond 270 Mpc we add contributions derived under the matter tracer model and the isotropy model (in doing so we assume that CR sources are numerous and distributed uniformly beyond 270 Mpc). In this procedure we first determine the fraction of the expected flux that comes from within 270 Mpc from eq. (9) (see also Fig. 5, top panel). We then randomly generate an event following either the isotropy model (with probability ) or the matter tracer model (with probability ).
In Sect. 4 we will be interested in the effect of source density on the statistical power of the test proposed in this work. In the case of a CR source density smaller than the density of galaxies in the KKKST catalog (i.e. ), we generate a reduced catalog with the desired source density by randomly selecting a fraction of galaxies from the full catalog. We then generate events from this reduced catalog following the procedure described above.
3 Anisotropy tests
3.1 The contrast parameter
The anisotropy test proposed in this work is based on the distribution of a statistical parameter within a set of CR events. We choose this parameter, which we denote as , to be equal to the expected flux from sources within the catalog distance in the matter tracer model, modulated with detector exposure:
where denotes the relative detector exposure, is defined in eq. (10), and are the (galactic) coordinates of the CR event, and is the threshold energy. Note that a choice of galaxy catalog is implicit in eq. (15) via eq. (10); we use the KKKST catalog unless stated otherwise. Because is a measure of the matter over- or underdensity in a given direction, we will refer to it as the contrast parameter. For the statistical tests considered in this work the normalization of is arbitrary.
In figure 7 we show the distribution of the contrast parameter when CR events are distributed isotropically and in the matter tracer model for three different CR energy thresholds. The figure clearly shows that takes preferentially larger values when events follow structure, as expected. This effect, which is most pronounced at the highest energies, underlies the statistical test proposed in the following section.
3.2 Test statistics
The statistical tests considered in this work are based on test statistics, real-valued numbers that characterize a set of CR events. In this section we introduce the four test statistics considered in this study. We will discuss the methodology of hypothesis testing in the following subsection.
First we consider the standard one-dimensional KS test statistic applied to the contrast parameter :
where denotes the cumulative distribution of in the test data set and denotes a reference cumulative distribution (all cumulative distributions will be assumed to be normalized). This reference distribution may be derived in the model that is to be tested, but it may also be derived in an alternative model.
Second, for the purpose of this study we propose a novel test statistic :
which represents a small modification of the standard KS test statistic . This statistic is also similar to the Kuiper test statistic, which is the sum of the two terms in eq. (17).
In order to compare the sensitivity of the - and -tests proposed here we also consider two test statistics that were discussed in the literature before. The first of these is the two-dimensional KS test applied to the galactic coordinates and of the CR events. We use the implementation of this test proposed in Ref. , which is based on earlier work in Ref. . The procedure put forward in Ref.  yields a quantity that approximates the largest difference between the reference and test distributions in the plane and over the four independent ways data can be accumulated in two dimensions:333Note that there are only two ways of accumulating data in the one-dimensional case, which are equivalent.
where denotes one of the four quadrants, with respect to a sampling point, in which data can be accumulated. Eq. (18) becomes an equality in the limit of an infinite number of sampling points. Due to computational constraints, we only consider a limited number of sampling points, namely those points contained in the test data set.
Finally we will also consider the test statistic that was recently proposed by Kashti and Waxman in Ref. , who found it to be more sensitive than methods based on the two-point correlation function or the angular power spectrum. We will denote this statistic as ( in Ref. ). The statistic is obtained after binning the sky and comparing the number of events per bin in the test data set to two independent model sets:
where denotes the observed number of events in bin in the test data; and denote the average numbers of events in bin in a reference model set and an alternative model set, respectively; and enumerates the bins.
3.3 Statistical power, significance, and -values
In this section we recall the definition and interpretation of statistical power, significance, and -values.
Our aim will be to test, i.e., to rule out or retain, a particular null hypothesis (e.g., “CR events are distributed isotropically”). Following the traditional frequentist approach put forward by Neyman & Pearson , one may distinguish two types of errors in this judgment:
|Type I error:||we reject while it is true;|
|Type II error:||we do not reject while an alternative hypothesis is true.|
The probabilities for the two types of error to occur are denoted as and , respectively. The first probability, , is known as the statistical significance. In a frequentist sense it is just the fraction of false rejections in many tests. Note that bears no reference to alternative hypotheses. The statistical power, defined as , is a measure of the ability of a given test to discriminate between the null hypothesis and an alternative hypothesis; it is computed for a particular alternative hypothesis.
The statistical significance and power are readily quantified for hypothesis tests that are based on the value of a single test statistic . The tests considered in this work fall in this category. In these tests one rules out the null hypothesis when falls in some critical range, which has been specified a priori. The significance is then equal to the probability that falls inside the critical range under ; the power of the test is the probability that falls inside the critical range when the alternative hypothesis is true. The significance can be determined from the distribution of under , denoted as , which may be known analytically or sampled by Monte Carlo event generation. Likewise, the power of the test is obtained from , the distribution of under . As an example consider the case that large values of are incompatible with . One would then decide to rule out when is in the critical range where is a predetermined critical value. The significance of this test is then given by the integral of over the interval and the power by the integral of over . In practice one generally determines from a certain required significance .
Different from statistical power or significance , the -value is a data-dependent quantity. It gives a figure of merit for the compatibility of a given measurement and the null hypothesis. More precisely, it is the probability that, under , the outcome of an experiment is at least as extreme as the observed outcome. Although small -values indicate greater incompatibility with , the -value has no precise frequentist interpretation. In the case of a hypothesis test based on a single test statistic, the -value of a given measurement can be determined by integrating over the range of -values that are equally or less compatible with the null hypothesis:
In the exemplary case that large values of correspond to poor agreement with the data, , where is the value of the test statistic corresponding to the observed data.
For the - and -test statistics introduced in the previous section, there is an important subtlety in choosing the parameter range that covers values of the test statistic “at least as extreme” as a particular value . This is because for these test statistics both very small and very large values indicate disagreement between the data and . As there is no a priori reason to reject only for very high or for very low values of the test statistic, we define as a combination of two disconnected parts that cover both very small and very large values of (two-sided interval). We choose these regions symmetrically such that they contain equal fractions of the test statistic distribution. This procedure is equivalent to defining the -value as follows:
where by definition when is higher than the average value and when is smaller.
An insightful way to relate -values to -probabilities is as follows. Equations (20) and (21) define a map . Therefore we can represent a measurement by its -value rather than its test statistic . Under the distribution of the -value (over many experiments) is, by definition, uniform between 0 and 1. We can thus rule out the null hypothesis, with significance , when we determine experimentally that , where is the -value associated with a measurement via eq. (20) and is a predetermined parameter. As a note of warning we would like to stress that it is inconsistent to interpret itself as a statistical significance, which amounts to adjusting to a posteriori.
4 Analysis of statistical power
In this section we consider the statistical power to rule out the “isotropy” model (CR events are distributed isotropically) when the matter tracer model holds. This power is a measure of the ability of the test to discriminate between the two models. We reiterate that, within the matter tracer model, we reconstruct the local matter density in the Universe from the KKKST catalog unless stated otherwise. We will, in an obvious notation, refer to the four hypothesis tests discussed in Sect. 3.2 as the -, -, -, and -tests. The results in this section are based on simulated event sets; we consider AGASA and PAO data in the next section.
We compute the statistical power to reject an isotropic distribution when the matter tracer model is true from the distributions of the test statistics under both models (see section 3.3). We sample these distributions by Monte Carlo event generation (see section 2.5). In this process we generate independent test event sets that each contain events. For every simulated event we compute the contrast parameter from eq. (15). We then compute the values of the test statistics from eqs. (16)(19) for every simulated test set. In these equations , and are computed from an isotropic reference event set, while is obtained from an alternative model event set generated in the matter tracer model. These reference and alternative model event sets contain () events. Because of the stochastic nature of our routine (Monte Carlo event generation), there is some numerical uncertainty in our results. We have verified that the results presented in this section are accurate within a few percent.
If the assumed CR source density is smaller than the density of galaxies in the KKKST catalog, we use different source realizations (reduced catalogs) to generate the test event sets under the assumption that CR sources trace the distribution of matter. Note that we compute the contrast parameter from the full catalog also in this case because in practice we would not know which fraction of the galaxies in the catalog are the actual CR sources.
We are now in the position to compute the statistical power of the hypothesis tests under consideration and investigate the dependence on the different parameters. To estimate the importance of event number and threshold energy we will consider three scenarios:
|100 CR events with energies above 40 EeV;|
|19 CR events with energies above 57 EeV|
Table 1 gives the statistical power of the four hypothesis tests under consideration to rule out an isotropic distribution when the matter tracer model is correct. These results correspond to the following values for the model parameters: source density fraction , injection spectral index , maximum injection energy eV (the exact value is unimportant in this regime), and smearing angle . Furthermore we have assumed no source evolution, and adopted the PAO detector exposure .
As may be seen in this table, the -statistic proposed in this work systematically yields a somewhat higher power than the standard KS test statistic . The power of the -test is much smaller than the power of the other tests in all three scenarios. Comparing the powers of the - and -tests, we find that the -test is superior in scenario I (large number of events at low energy) while the -test performs better in scenario III (small number of events at high energies). In the intermediate cases II (PAO data) the powers are comparable, with the -test performing slightly better.
We now focus on the -test, the most sensitive of the two tests proposed in this paper, and investigate how the statistical power changes when we vary model parameters. This is indicated in tables 2 and 3. The second column of these tables, labeled “reference”, corresponds to the reference values for the parameters which were also used in table 1. The other columns indicate how the statistical power changes when we vary one of these parameters, and when we use the PSCz catalog rather than the KKKST catalog.
First let us consider the effect of the CR source density. As table 2 indicates, the statistical power of the test decreases with decreasing source density. This is due to the fact that the distribution of when sources trace the distribution of matter becomes broader (being averaged over many source realizations), thus increasing the overlap with the distribution of when CRs are distributed isotropically. This is illustrated in figure 8 for scenario II. As may be seen from the table, the statistical power convergences to an asymptotic value as approaches 1. This confirms our earlier estimates that the many-source regime applies for a CR source density equal to the density of galaxies in the KKKST catalog.
In table 3 we show how the statistical power varies when we change different model parameters ( here). First, we have varied the power-law index for the injected proton spectrum. We observe that the statistical power decreases when the injection spectrum becomes harder. This may be understood from the fact that the effective viewing distance increases with a harder injection spectrum. The increase in observable volume dilutes the contrast from nearby structure. Second, we consider the effect of maximum injection energy. Here we find that the statistical power increases when the maximum CR energy upon injection is not much larger than the threshold energy. This can also be understood in terms of viewing distance, with a lower implying a smaller viewing distance and hence a stronger contrast. Third, we have increased the smearing angle to . This obviously dilutes anisotropies and hence decreases the statistical power to rule out an isotropic distribution when the matter tracer model is true. Fourth, we have investigated the statistical power of the proposed test when the distribution of matter is derived from the PSCz catalog  instead of the KKKST catalog. Here we weight the flux of individual galaxies in the catalog with the inverse of the selection function (see Ref. ), adjusted to km s Mpc. To avoid large fluctuations due to a few far-away sources, we cut the PSCz catalog at 270 Mpc and assume an isotropic distribution at larger distance. This distance is chosen for comparison with the KKKST catalog. For the same reason we also cut the PSCz catalog at . We observe that the statistical power decreases substantially when we replace the KKKST catalog by the PSCz catalog. This is due to the fact that matter density contrasts are less pronounced in the PSCz catalog, an issue which will be discussed in Sect. 6.2. Finally (not shown in the table), we find that source evolution has little effect on the statistical power. Even when we assume that CR sources trace strong AGN luminosity density evolution as parameterized in Ref. , the results remain the same within numerical uncertainties. This is expected because high-energy CRs originate from sources in the local Universe, rendering the effect of source evolution small.
We have compared the statistical powers of the - and -test in the first part of this section. Qualitatively, we have found that the -test performs better at high energies and the -test better at low energies. The precise energy where the two tests are equally powerful depends on many parameters, such as the experiment’s field-of-view, the deflection angle, and the assumed CR spectrum. Using reference values of the model parameters we found that the -test is the most sensitive in scenario III, while the -test is the most sensitive for scenarios I and II (see table 1). We have verified explicitly that this conclusion remains valid for the different parameter values that were considered in tables 2 and 3.
4.3 Comparison with results in the literature
Refs. [9, 11, 13] explicitly give the statistical power to rule out an isotropic CR distribution when CR sources follow structure for other statistical tests, thus allowing for a comparison with our estimates. In table 4 we compare the statistical power of the -test proposed here to the results of Refs. [9, 11, 13]. In obtaining our estimates we vary event numbers, threshold energies, statistical significance, and other model parameters to match the values that were used in the quoted references. (We however stick to the KKKST catalog for the distribution of galaxies). Results between the different references should therefore not be compared. The results in the table show that the power of the -test, as determined in this work, is comparable to or larger than the power of the tests proposed in the quoted references.
|60||31||PAO||0.05||0.56||0.25||(; Mpc case)|
The power estimates in table 4 depend both on the test itself and on the modeling of CR arrival directions. If one is interested in a comparison of the tests per se, the tests should be applied to exactly the same event sets. We have done this for the -statistic put forward in Ref.  (denoted as in this paper) in the previous subsections; table 1 shows how this test statistic compares to the -test in performance.444We note that the power of the -test for 100 events above 40 EeV as reported in table 1 () is significantly larger than the results given in Ref.  (). This is primarily due to the different catalogs that were used; we will come back to this in Sect. 6.2.
5 Anisotropy in AGASA and PAO data
In this section we assess the compatibility of the publicly available data from AGASA [41, 42] and PAO  with an isotropic distribution and predictions of the matter tracer model, with the local matter density in the Universe reconstructed from the KKKST catalog. The data sets are cut at , resulting in a total of 44 events above 40 EeV for AGASA and a total of 19 events above 57 EeV for PAO. The results presented here are primarily intended to illustrate the method presented in this work, the number of available events still being small and the energy calibration uncertain.
|value of test statistic in AGASA data||0.087||0.050||-0.59||0.19|
|-value for “isotropy” model||0.86||0.61||0.88||0.20|
|-value for matter tracer model||0.98||0.20||0.018||0.45|
|value of test statistic in PAO data||0.25||0.18||6.3||0.16|
|-value for “isotropy” model||0.15||0.20||0.035||0.83|
|-value for matter tracer model||0.58||0.66||0.85||0.88|
The -values corresponding to the four statistical tests under consideration are determined by integrating the test statistic distributions in the model to be tested over the range of possible outcomes that are at least as extreme as the observed one (see section 3.3). The test statistics are determined from eqs. (16)(19), where we compute the quantities , , and in the isotropy model and the quantity in the matter tracer model. (This is a matter of choice, and one could also consider a test where the reference quantities are computed in the matter tracer model. We do not pursue this possibility here.) The numerical sampling of the test statistic distributions is discussed in Sect. 4.1.
The exposure of both AGASA and PAO is essentially geometrical at the highest energies. We model this exposure following Ref. , adopting a detector latitude () and a maximum viewing angle () for AGASA (PAO).
In table 5 we show the values of the -, -, - and -test statistics for the AGASA data [41, 42], together with the associated -values for the “isotropy” model and the matter tracer model. The -values are all above 0.01 and are therefore not conclusive. The general trend is that the data are more compatible with isotropy than with the matter tracer model. In particular, the -test statistic yields a -value of for the AGASA data and the matter tracer model, indicating mild incompatibility. These results are compatible with the absence of large-scale anisotropies reported in Ref. .
Table 6 shows the values of the test statistics for the PAO data  together with the -values for the “isotropy” and matter tracer model. Also here the -values are larger than 0.01 and hence inconclusive. All tests under consideration are fully compatible with the matter tracer model. The -test shows a very mild () incompatibility between the PAO data and the isotropy model predictions, while the other tests are well compatible with isotropy.
It is not very surprising that an isotropic distribution cannot be ruled out with a strong significance with 19 events above 57 EeV when the matter tracer model is true: we have estimated in the previous section that the probability to find a -value smaller than 0.01 for the PAO data and the isotropy model ranges from 3.3% (-test) to 26% (-test) in this case (see table 1).
We expect that the tests proposed in this work will constrain models of CR origin when more data become available. We estimate that events above 57 EeV are required to have a 50% probability of excluding an isotropic distribution with significance 0.01 if the matter tracer model is correct. (We obtained this results with the procedure discussed in the previous section). The required data can be accumulated by PAO and the Telescope Array (TA) experiment  within a few years. A determination of the optimal search strategy (including optimization of the energy threshold) requires an extensive study of the parameter space and is beyond the scope of this paper.
We have shown in Sect. 4 that the - and -tests are very efficient in discriminating between the matter tracer model and an isotropic distribution of UHECRs, especially for the highest energy part of the UHECR spectrum where the number of events is small and the flux map exhibits more contrast. The sensitivity comes at a price: one may imagine deviations from the matter tracer model to which our test is blind. The existence of such deviations is clear from Fig. 1 (which corresponds to the example discussed in the introduction): the histogram on the left does not change if one moves the events leaving them on the same flux contour (i.e. the same grayscale). More generally, the - and -test statistics as defined in eqs. (16) and (17) are unchanged if we move CR events along contours of equal . In particular, one may group all events that fall on the same flux contour in one place without changing the - and -test statistics. Although this would certainly increase the anisotropy, it would not be noticed by these tests. The -test proposed in Ref.  suffers from the same kind of blindness, as it is invariant under redistributions of events over bins with equal fraction .
6.2 Comparison with earlier results for the -test
Here we compare our results for the statistical power and -values of the -test with results presented by Kashti and Waxman in Ref. , where this statistic was introduced.555It is denoted as in Ref. . We compare our results with the largest source density considered in Ref. , because we are interested in the many-source regime, and with the “unbiased” model because this resembles our setup most closely.
Our estimates for the statistical power to reject an isotropic distribution when the matter tracer model is true are significantly larger than the estimates in Ref. , which is to say that we find a much better discrimination between the two models. For scenario I ( events above EeV) we find a power (significance 5%; see table 1), whereas Kashti and Waxman find (same significance; their table II). We have verified that our estimates are also larger for other values of and for which power estimates are given in Ref. ; these do not include our scenarios II and III though. Although our model assumptions differ in more than one respect, the observed difference appears to be primarily due to the different catalogs that were used in modeling CR arrival directions, namely the PSCz catalog  in Ref.  and the KKKST catalog in this work (Sect. 2.1). To show this we have recomputed the statistical power for the -test using the PSCz catalog and find that the power reduces to . When we also change model parameters and cuts to match the values used in Ref.  (i.e., , strong source evolution, cut the galactic plane above ) we find , in reasonable agreement with the result obtained in Ref. . (Comparing our results to data files referred to in Ref. , we trace the fact that our estimates are somewhat higher to a stronger contrast in the matter densities derived from the PSCz catalog). The relatively strong impact of the source catalog on the statistical power is in keeping with results obtained for the -test in section 4 (see table 3).
Now we compare the -values for an isotropic model distribution and the PAO data between this study and Ref. . Our -value of (see table 6) is almost an order of magnitude larger than the value of reported in Ref.  (their table III), meaning that we find a much better agreement between the data and an isotropic CR distribution. At first sight this may appear to be in contradiction with the fact that our power estimates are larger than those in Ref. : the larger statistical power leads us to expect — if the matter tracer model describes nature correctly — that we would find a stronger incompatibility between “isotropy” and PAO data (more precisely, we have better chances of finding a strong incompatibility). This discrepancy can be explained if the true CR source density contrasts are larger than the density contrasts in the PSCz catalog used in Ref. , so that the model predictions based on that catalog are not accurate. This may (partly) be due to selection effects, as it has been pointed out that the limited angular resolution of the IRAS instrument used to compile the catalog may lead to a systematic underestimate of the number density of galaxies in high-density regions (see discussion in Ref. ). It may also indicate that CR sources are intrinsically biased with respect to the infrared galaxy sample underlying the PSCz catalog. The results presented in Ref.  lend further support to the hypothesis that CR source density contrasts are larger than those derived from the PSCz catalog: The -value for the PAO data and the “unbiased” model (in which CR sources trace the distribution of matter) in Ref.  is also small, viz. (table III in their work), which indicates a mild discrepancy between the PAO data and the model predictions based on the PSCz catalog. (We find a -value of 0.85, and hence perfect agreement, when we apply the -test to the PAO data and the matter tracer model using the KKKST catalog; see table 6). In fact, the PAO value of the -test statistic in Ref.  is larger than the average value predicted in their “unbiased” model, which in turn is larger than the average value predicted under the assumption of an isotropic event distribution. This suggests that the data exhibit more contrast than the “unbiased” model predictions. From the same observation, the authors of Ref.  conclude that the PAO data show a preference towards their “biased” model, which has more contrast than the “unbiased” model. Future CR data will provide further clues as to the actual distribution of sources and the relative strength of over- and underdensities.
In this paper we have presented a new and powerful method to test models of cosmic-ray origin that predict a continuous distribution of CRs over the sky. A prominent model of this kind is the “matter tracer model”, a simple yet general model that predicts that local matter overdensities are reflected in the distribution of CRs. Its main assumptions, a continuous distribution of sources tracing the distribution of matter in the Universe and small CR deflections, represent the limiting case of many detailed scenarios. We therefore consider this model as a useful benchmark model in CR anisotropy studies. The local matter density in the Universe, which is required as input for the matter tracer model, is derived from a galaxy catalog. For this purpose we primarily use the catalog compiled by Kalashev et al. , a complete sample of galaxies up to 270 Mpc (the KKKST catalog; see Sect. 2.1). We have however recomputed some of our results using the PSCz catalog , which was used in various earlier studies (e.g., Ref. ), for comparison.
We analyzed in detail (Sect. 4) the statistical power of the two proposed tests (denoted as the - and -tests) to reject an isotropic CR distribution when the matter tracer model is true (with matter densities reconstructed from the KKKST catalog). We find that the proposed tests are competitive or superior in power to other tests in the literature. They are especially sensitive in the physically interesting case of relatively few events with very high energies, viz. events above EeV (the test statistic proposed in Ref.  provides a more powerful discrimination in the case of many events with lower threshold energy, EeV). The proposed tests are completely binless, which is an important advantage because binning introduces an ambiguity that complicates the statistical analysis.
We have applied our tests to the AGASA data above 40 EeV [41, 42] and the PAO data above 57 EeV , and found that the present data are compatible with both the matter tracer model (with matter densities reconstructed from the KKKST catalog) and with an isotropic CR distribution (Sect. 5). The general trend is that the AGASA data show a greater compatibility with an isotropic distribution, whereas the PAO data favor the matter tracer model. We expect that the proposed tests will provide meaningful constraints on models when more data become available. We estimate that events above 57 EeV are required to have a 50% probability of excluding an isotropic distribution with significance 0.01 if the matter tracer model is correct. The required data can be accumulated by PAO and TA within a few years.
Summarizing our method, both the - and - tests are based on the (continuous) distribution of a parameter — which is a measure of the amount of matter in the nearby Universe — over a set of observed CRs (Sect. 3.1). In the matter tracer model the probability that a CR comes from an overdense region is higher than from an underdense region, so that the -parameter takes preferentially high values. For an isotropic distribution there is no such correlation and is essentially a random variable. Correlating the distribution of in a set of observed CR events with model distributions then gives a measure of the (in)compatibility between observations and model predictions. This last step may be done with the standard Kolmogorov-Smirnov test (-test) or a modified version thereof proposed in this work (-test; Sect. 3.2). The latter is shown to be somewhat more powerful in discriminating between the matter tracer model and an isotropic CR distribution.
Reconstructing the matter density in the Universe from a galaxy catalog, as done in this study, necessarily introduces some ambiguity and uncertainties. For the KKKST catalog, which we have primarily used, the dominant source of errors is in the photometric redshift estimates for galaxies beyond 30 Mpc. The associated uncertainties are acceptable given the other uncertainties we are faced with, most notably in the energy resolution of present CR experiments and in modeling CR energy loss in the Universe. For the purpose of this study, the most important difference between the KKKST catalog and the PSCz catalog is that the former exhibits stronger density contrasts, which improves the prospects of discriminating between the matter tracer model and an isotropic CR distribution. Whether or not actual CR sources exhibit such a strong density contrast remains an open question, to be answered by observations. We would like to point out that the increased contrasts cannot be an artefact of the random errors associated with photometric redshift estimates. Such errors will, on average, bring sources closer. Keeping the CR horizon fixed, this results in a more homogeneous flux distribution rather than a more contrasted one (as may be easily seen by noting the equivalence with increasing the horizon while keeping distances fixed).
Current and future data can test model predictions based on a given galaxy catalog. Using the present data from PAO, we argued in Sect. 6.2 that the CR source density seems to exhibit stronger contrasts than what is derived from the PSCz catalog, while the data are fully compatible with the KKKST catalog. This picture is however not supported by the AGASA data, which favor an isotropic model. In the future it may become possible to investigate the spatial distribution of CR sources in a greater amount of detail by comparing CR data with predictions from catalogs based on various galaxy types with different clustering properties. This would require more, and more precise, CR data as well as improved galaxy catalogs.
We comment on some possible improvements to and extensions of the proposed tests. A first possibility is to include the energies of individual events rather than a single threshold energy in the analysis. A preliminary investigation has shown that this can improve the statistical power, but more work is needed to study this idea in detail. Second, we note that it is possible to choose the statistical parameter in a different manner than eq. (15), e.g. to be equal to the detector exposure (so that will take preferentially high values in an isotropic model). We have not investigated this possibility here. Finally the estimates for proton energy-loss lengths may be updated using more sophisticated numerical propagation routines when future data call for a higher degree of precision.
Acknowledgments.We thank Sergey Troitsky for valuable discussions and for providing us with the KKKST catalog, and Alessandro Cuoco, Michael Kachelrieß, Tamar Kashti, and Eli Waxman for comments. H.K. would like to thank Thomas Lessinnes for illuminating discussions on hypothesis testing. H.K. and P.T. are supported by Belgian Science Policy under IUAP VI/11 and by IISN. The work of P.T. is supported in part by the FNRS, contract 1.5.335.08.
-  HiRes Collaboration, R. Abbasi et. al., Observation of the GZK cutoff by the HiRes experiment, Phys. Rev. Lett. 100 (2008) 101101, [astro-ph/0703099].
-  Pierre Auger Collaboration, J. Abraham et. al., Observation of the suppression of the flux of cosmic rays above eV, Phys. Rev. Lett. 101 (2008) 061101, [arXiv:0806.4302].
-  Pierre Auger Collaboration, J. Abraham et. al., Correlation of the highest energy cosmic rays with nearby extragalactic objects, Science 318 (2007) 938, [arXiv:0711.2256].
-  Pierre Auger Collaboration, J. Abraham et. al., Correlation of the highest-energy cosmic rays with the positions of nearby active galactic nuclei, Astropart. Phys. 29 (2008) 188, [arXiv:0712.2843].
-  K. Dolag, D. Grasso, V. Springel, and I. Tkachev, Constrained simulations of the magnetic field in the local universe and the propagation of UHECRs, JCAP 0501 (2005) 009, [astro-ph/0410419].
-  G. Sigl, F. Miniati, and T. A. Ensslin, Ultra-high energy cosmic ray probes of large scale structure and magnetic fields, Phys. Rev. D70 (2004) 043007, [astro-ph/0401084].
-  S. L. Dubovsky, P. G. Tinyakov, and I. I. Tkachev, Statistics of clustering of ultra-high energy cosmic rays and the number of their sources, Phys. Rev. Lett. 85 (2000) 1154, [astro-ph/0001317].
-  M. Kachelriess and D. Semikoz, Ultra-high energy cosmic rays from a finite number of point sources, Astropart. Phys. 23 (2005) 486, [astro-ph/0405258].
-  T. Kashti and E. Waxman, Searching for a Correlation Between Cosmic-Ray Sources Above eV and Large-Scale Structure, JCAP 0805 (2008) 006, [arXiv:0801.4516].
-  E. Waxman, K. B. Fisher, and T. Piran, The signature of a correlation between -eV cosmic ray sources and large scale structure, Astrophys. J. 483 (1997) 1, [astro-ph/9604005].
-  A. Cuoco, R. D. Abrusco, G. Longo, G. Miele, and P. D. Serpico, The footprint of large scale cosmic structure on the ultra-high energy cosmic ray distribution, JCAP 0601 (2006) 009, [astro-ph/0510765].
-  M. Kachelriess and D. V. Semikoz, Clustering of ultra-high energy cosmic ray arrival directions on medium scales, Astropart. Phys. 26 (2006) 10, [astro-ph/0512498].
-  O. E. Kalashev, B. A. Khrenov, P. Klimov, S. Sharakin, and S. V. Troitsky, Global anisotropy of arrival directions of ultra-high- energy cosmic rays: capabilities of space-based detectors, JCAP 0803 (2008) 003, [arXiv:0710.1382].
-  P. Sommers, Cosmic Ray Anisotropy Analysis with a Full-Sky Observatory, Astropart. Phys. 14 (2001) 271, [astro-ph/0004016].
-  S. Singh, C.-P. Ma, and J. Arons, Gamma-ray bursts and magnetars as possible sources of ultra high energy cosmic rays: Correlation of cosmic ray event positions with IRAS galaxies, Phys. Rev. D69 (2004) 063003, [astro-ph/0308257].
-  A. Smialkowski, M. Giller, and W. Michalak, Luminous infrared galaxies as possible sources of the UHE cosmic rays, J. Phys. G28 (2002) 1359, [astro-ph/0203337].
-  N. W. Evans, F. Ferrer, and S. Sarkar, The anisotropy of the ultra-high energy cosmic rays, Astropart. Phys. 17 (2002) 319, [astro-ph/0103085].
-  D. Harari, S. Mollerach, and E. Roulet, Kolmogorov-Smirnov test as a tool to study the distribution of ultra-high energy cosmic ray sources, Mon. Not. Roy. Astron. Soc. 394 (2009) 916, [arXiv:0811.0008].
-  H. Takami, T. Nishimichi, K. Yahata, and K. Sato, Cross-Correlation between UHECR Arrival Distribution and Large-Scale Structure, arXiv:0812.0424.
-  W. Saunders et. al., The PSCz Catalogue, Mon. Not. Roy. Astron. Soc. 317 (2000) 55, [astro-ph/0001117].
-  T. H. Jarrett et. al., 2MASS Extended Source Catalog: Overview and Algorithms, Astron. J. 119 (2000) 2498, [astro-ph/0004318].
-  G. Paturel, C. Petit, P. Prugniel, G. Theureau, J. Rousseau, M. Brouty, P. Dubois, and L. Cambrésy, HYPERLEDA. I. Identification and designation of galaxies, Astron. Astrophys. 412 (2003) 45.
-  T. Jarrett, Large Scale Structure in the Local Universe: The 2MASS Galaxy Catalog, astro-ph/0405069.
-  D. Fargion, Light Nuclei solving Auger puzzles?, Phys. Scripta 78 (2008) 045901, [arXiv:0801.0227].
-  S. Gureev and S. Troitsky, Physical conditions in potential sources of ultra-high- energy cosmic rays. II. Nearby active galaxies correlated with Auger events, arXiv:0808.0481.
-  M. J. Hardcastle, C. C. Cheung, I. J. Feain, and L. Stawarz, High-energy Particle Acceleration and Production of Ultra- high-energy Cosmic Rays in the Giant Lobes of Centaurus A, Mon. Not. Roy. Astron. Soc. 393 (2009) 1041, [arXiv:0808.1593].
-  D. Gorbunov, P. Tinyakov, I. Tkachev, and S. V. Troitsky, Comment on ‘Correlation of the Highest-Energy Cosmic Rays with Nearby Extragalactic Objects’, JETP Lett. 87 (2008) 461, [arXiv:0711.4060].
-  T. Wibig and A. W. Wolfendale, Heavy Cosmic Ray Nuclei from Extragalactic Sources above ‘The Ankle’, arXiv:0712.3403.
-  I. V. Moskalenko, L. Stawarz, T. A. Porter, and C. C. Cheung, On the Possible Association of Ultra High Energy Cosmic Rays with Nearby Active Galaxies, Astrophys. J. 693 (2009) 1261, [arXiv:0805.1260].
-  Pierre Auger Collaboration, M. Unger, Study of the Cosmic Ray Composition above 0.4 EeV using the Longitudinal Profiles of Showers observed at the Pierre Auger Observatory, arXiv:0706.1495. Contribution to the 30th ICRC.
-  A. V. Glushkov et. al., Muon content of ultra-high-energy air showers: Yakutsk data versus simulations, JETP Lett. 87 (2008) 190, [arXiv:0710.5508].
-  The High Resolution Fly’s Eye Collaboration, R. U. Abbasi et. al., A Study of the Composition of Ultra High Energy Cosmic Rays Using the High Resolution Fly’s Eye, Astrophys. J. 622 (2005) 910, [astro-ph/0407622].
-  V. Berezinsky, A. Z. Gazizov, and S. I. Grigorieva, On astrophysical solution to ultra high energy cosmic rays, Phys. Rev. D74 (2006) 043005, [hep-ph/0204357].
-  H. B. J. Koers and P. Tinyakov, Relation between the neutrino flux from Centaurus A and the associated diffuse neutrino flux, Phys. Rev. D78 (2008) 083009, [arXiv:0802.2403].
-  M. Kachelriess, E. Parizot, and D. V. Semikoz, Constraining the cosmic ray source spectrum from observations in the GZK regime, arXiv:0711.3635.
-  C. D. Dermer, S. Razzaque, J. D. Finke, and A. Atoyan, Ultra High Energy Cosmic Rays from Black Hole Jets of Radio Galaxies, arXiv:0811.1160.
-  G. Fasano and A. Franceschini, A multidimensional version of the Kolmogorov-Smirnov test, Mon. Not. Roy. Astron. Soc. 225 (1987) 155.
-  J. A. Peacock, Two-dimensional goodness-of-fit testing in astronomy, Mon. Not. Roy. Astron. Soc. 202 (1983) 615.
-  J. Neyman and E. S. Pearson, On the problem of the most efficient tests of statistical hypotheses, Phil. Trans. R. Soc. A 231 (1933) 289.
-  B. J. Boyle and R. J. Terlevich, The cosmological evolution of the QSO luminosity density and of the star formation rate, Mon. Not. Roy. Astron. Soc. 293 (1998) L49, [astro-ph/9710134].
-  M. Takeda et. al., Small-scale anisotropy of cosmic rays above -eV observed with the Akeno Giant Air Shower Array, Astrophys. J. 522 (1999) 225, [astro-ph/9902239].
-  N. Hayashida et. al., Updated AGASA event list above eV, astro-ph/0008102.
-  H. Tokuno et. al., The Telescope Array experiment: Status and prospects, J. Phys. Conf. Ser. 120 (2008) 062027.