© ESO, 2014

1. Introduction

Discovering the origin of the ultra-high-energy cosmic rays (UHECRs), with energies on the order of 100 EeV is widely recognized as an important challenge in the field of astroparticle physics, both from the theoretical and observational points of view. Among the main scientific goals are a deeper understanding of particle acceleration in the universe and the search for complementary information and constraints on potential sources, in a multimessenger approach involving photons, energetic nuclei, neutrinos, and perhaps gravitational waves (see Kotera & Olinto 2011, for a recent review).

This quest has been pursued for decades with larger and larger detectors (see Nagano & Watson 2000; Letessier-Selvon & Stanev 2011, for reviews) and is currently led in the southern hemisphere by the Pierre Auger Collaboration (Pierre Auger Collaboration 2004) and in the northern hemisphere by the Telescope Array Collaboration (Abu-Zayyad et al. 2012a; Tokuno et al. 2012). Important milestones have been reached, notably through the confirmation of the so-called GZK effect (Greisen 1966; Zatsepin & Kuz’min 1966), which causes a rapid decrease in the UHECR flux around 60 EeV (Abbasi et al. 2008; Pierre Auger Collaboration 2008b, 2010b; Abu-Zayyad et al. 2013a) owing to the interaction of the UHE protons and/or nuclei with the extragalactic background radiation. Hints of anisotropies in the arrival directions of the UHECRs above ~60 EeV have also been reported (Pierre Auger Collaboration 2007, 2008a), indicating that the intervening magnetic fields do not completely randomize the distribution of UHECRs on the sky, as it does at lower energies, preventing a direct identification of the sources.

However, the initial hope of revealing the sources by observing an accumulation of UHECR events in well defined directions has been frustrated up to now. This is either because the number of contributing sources is too large and very few events have been observed from any given source, or because the particles observed from a given source are deflected and spread over areas of the sky that are too large, or both.

To overcome these difficulties, a natural strategy is to focus on the highest energy cosmic rays, because their magnetic rigidity is expected to be greater (unless they also have a higher charge), and because the number of sources visible from Earth is much smaller, due to the GZK effect. General statistical studies taking the relevant propagation effects into account show that only a handful of sources are responsible for most of the detectable UHECRs around 100 EeV. Depending on the source density and the exact contingent distribution of the sources around us, the dominant source in the sky typically accounts for between 15% and 60% of the total flux at 100 EeV (Blaksley et al. 2013). In such conditions, it seems likely that the very few contributing sources could be distinguished as separate clusters of events on the sky, even if most particles suffer relatively large deflections, provided enough UHECRs are observed.

This has not yet occurred with the very limited statistics available. No significant anisotropy has been detected at 100 EeV, and even though the Pierre Auger Observatory has reported a departure from isotropy at lower energy, E ~ 60 EeV, at the 99% confidence level (Pierre Auger Collaboration 2007, 2010c), its interpretation is still unclear, and it could not be used to gather any decisive information about the sources.

In total, an exposure ~30 000 km² sr yr has been accumulated so far, and it appears that decisive progress toward identifying the sources will not be possible unless a significant increase in exposure is achieved. In this paper, we investigate whether a next-generation detector gathering an exposure of 300 000 km² sr yr at the highest energies can indeed detect a significant anisotropy in the UHECR arrival directions, even though no clear signal emerges with the currently available statistics.

To this end, we investigate several astrophysical scenarios, varying the UHECR source density as well as the injection composition and energy spectrum, and build simulated sky maps taking into account the propagation of the UHECRs in the cosmological photon backgrounds and in the extragalactic and Galactic magnetic fields. These sky maps, built for various total numbers of events, are then analysed from the point of view of their intrinsic anisotropy, through the analysis of the two-point correlation function. Among the different scenarios, only those that lead to typical sky maps that are not ruled out by the current data (i.e., when examined with the Pierre Auger Observatory statistics) are considered as reasonable scenarios, and explored with more statistics. Many realizations of each astrophysical scenario are simulated, to investigate the range of possible corresponding sky maps, referred to as the “cosmic variance”.

To be definite, we use the JEM-EUSO (the Extreme Universe Space Observatory on board the Japanese Experiment Module of the International Space Station) mission as a prospective example of such a next generation detector, i.e., we assume a (nearly) uniform full-sky coverage and the actual detection efficiency of JEM-EUSO as a function of energy, as described in Adams et al. (2013). We also assume a conservative energy resolution of 30% to demonstrate that the resulting spillover of lower energy events (whose sources can be more distant and numerous) due to incorrectly reconstructed energy does not compromise the main result of this study: significant anisotropy is indeed expected to be detected at the highest energies for essentially all the models studied, even in the extreme case where the UHECRs are completely dominated by heavy nuclei at the highest energies. Therefore, a full-sky coverage detector achieving an exposure of 300 000 km² sr yr should make a valuable difference in the current state of knowledge regarding UHECRs and start a new phase in their phenomenological and theoretical study by identifying significant anisotropies and studying individual sources.

In Sect. 2, we emphasize the key feature of the GZK horizon effect, which makes such a study possible at the highest end of the cosmic ray energy spectrum, even in the presence of relatively large deflections (e.g., with an Fe-dominated composition). In Sect. 3, we describe the various astrophysical scenarios explored and their main ingredients. In Sect. 4, we describe the magnetic field models implemented in the simulation and the general procedure used to produce representative UHECR sky maps. In Sect. 5, we present and discuss the main results of our investigation, through the quantitative study of the anisotropies that can be expected for the various scenarios. A general summary is given in Sect. 6.

2. The GZK cut-off and UHECR anisotropies

As recalled in the introduction, current UHECR observatories have not succeeded in detecting significant and unambiguous anisotropies or in reaching the stage of individual source astronomy. Two important elements of the UHECR phenomenology may concur in such a situation: the typical source density and the angular deflection of the particles. Since the total UHECR flux is known, the effective source density is directly related to the intrinsic power of the sources, which determines how many events can be expected from a given source (depending on its distance). As Blaksley et al. (2013) have shown, the most intense sources should have already contributed several events to the current dataset. Nevertheless, no obvious multiplets, defined as UHECR events coming from the same source, have been identified so far. If multiplets are present, but not recognized as such, it must be that the typical deflections exceed or approach the angular separation between sources. At energies between 10 EeV and 50 EeV, say, this may be due to the large number of sources, which overlap more or less uniformly over the sky. At higher energies, however, the GZK effect can be extremely valuable because of its most basic consequence, which is to reduce the horizon distance; i.e., the distance beyond which the contribution of UHECR sources to the observed flux is negligible. At the same time as the horizon becomes significantly closer with increasing energy (from ~180 Mpc to ~75 Mpc between 60 EeV and 100 EeV for protons), the number of visible sources is considerably reduced (by more than a factor of 10 between these two energies), and the background of more distant sources – which are also more isotropically distributed – is cut.

For this reason, it is crucial to focus on the highest energy particles. Above 60 EeV or so, the visible sources are located in a region of the universe where matter is not distributed uniformly. Even in the case of high source densities, the source distribution itself would leave a visible imprint on the UHECR sky map if the deflections were small (Decerprit et al. 2012). Thus, the absence of clear large-scale anisotropy patterns already puts some constraints on the typical deflection of UHECRs, with deflections probably greater than ~10–15 degrees, whether from a high magnetic field or a large electric charge. This agrees with the indication that a transition toward heavier nuclei occurs above 10¹⁹ eV, as shown by the analysis of the Pierre Auger Observatory data (Pierre Auger Collaboration 2010a). A small fraction of low-Z nuclei may nevertheless be present among the UHECRs, which would eventually lead to the observation of multiplets on a small angular scale as the statistics increases, especially at the highest energies where the particle deflections are smallest.

To overcome this general problem, a natural strategy is thus to extend the exposure of the experiments at energies two or three times higher than currently available with reasonable statistics, not only to reduce the deflections by the same factor (if the composition does not become heavier above 60 EeV), but mostly to drastically reduce the number of sources and make their separation on the sky easier. At 100 EeV, most of the UHECRs are expected to come from a handful of sources, even in the case of relatively high source densities (Blaksley et al. 2013). One may thus expect to be able to isolate them or at least detect a significant anisotropy from the associated multiplets, even if the deflections are relatively large, e.g., in the case of Fe nuclei in a typical Galactic and extragalactic magnetic field.

To verify this basic idea quantitatively and estimate the required integrated exposure, we perform simulations under various astrophysical assumptions, as described in the next section.

The reason many different scenarios are possible is that the all-sky UHECR spectrum does not contain much information by itself, with the current statistics. The significant, drastic reduction of the flux observed above ~60 EeV is generically attributed to “the” GZK cut-off. However, there is not a single GZK cut-off, and the current measurements are compatible with a wide range of models, because what gives its shape to the all-sky propagated spectrum is the horizon scale structure (its variation with energy), much more than the individual source spectrum, maximum energies, source density, and even source composition. The coincidental similarity, in first approximation, of the horizon scale structure for protons and Fe nuclei (Allard et al. 2008) is such that the current spectrum can be accounted for by a model that assumes pure proton sources, pure Fe sources, or a mixed composition of cosmic rays. Even unrealistic scenarios where the sources accelerate, say, only C nuclei, only Si nuclei, or essentially any individual nuclear species, can give a good fit of the data (above the ankle, interpreted as the transition from Galactic to extragalactic nuclei in these scenarios), provided that one suitably adjusts the source spectrum, which is unknown anyway. This is because Fe nuclei have the same horizon structure as the protons, and the lighter nuclei are rapidly destroyed by photodisintegration in the intergalactic photon fields. The resulting secondary protons then give rise to the standard GZK cut-off, as if the sources were effectively proton sources, i.e., as if they were essentially emitting the secondary particles directly.

This property of UHECR phenomenology is why we do not expect significant progress in this field of research from a refined study of the all-sky spectrum, but rather from anisotropy studies at the highest energies and/or from the comparison of the energy spectra in different regions of the sky. However, one should also keep in mind that the discovery of a recovery in the energy spectrum (see Berezinsky et al. 2006, for instance) or any unexpected feature would be significant.

To determine whether significant anisotropies can be observed with future detectors and to evaluate the confidence level with which they can be expected to be measured, we simulated UHECR sky maps for a wide range of models, as described below.

3. The source models and their parameters

3.1. General astrophysical assumptions

Standard propagation studies that take the energy losses and angular deflections of the UHECRs into account allow the whole sky spectrum to be reproduced with a limited number of parameters, under the simplifying assumption that all UHECR sources are essentially identical; i.e., i) they have the same intrinsic power; ii) they inject UHECRs into the extragalactic medium continuously; iii) they do so with the same power-law energy spectrum; and iv) with the same composition. The remaining free parameters are the source density, the logarithmic slope of the source spectrum, the maximum energy of the UHECRs at the source, and the relative abundances of the various nuclei. These parameters are not independent and must be chosen so as to reproduce the observed spectrum (see below).

It is likely that, in reality, individual sources are all different. However, the above assumptions allow a large set of models to be explored that should be representative of the range of patterns one may expect from the point of view of anisotropies. Relaxing them would introduce more free parameters on which there are no constraints at the moment, whether observationally or theoretically, without significantly enlarging the range of possibilities explored. The only assumption that we relax in the present study is that of an identical intrinsic power for all sources. We adopt the same intrinsic luminosity distribution as that of the galaxies in the catalog that we use to represent a realistic source distribution in the local universe (see below). This can be argued to be providing a more natural benchmark luminosity distribution than the “standard candle” assumption, but it is a non-essential feature of our study, so we stick to this initial assumption throughout and do not explore an extra dimension of the parameter space by varying the luminosity distribution.

A distribution of maximum energies can also be expected in principle, but its net effect would be a difference between the actual spectrum of the UHECRs at the source and the effective source spectrum convolving the former with the maximum energy distribution, in a way that does not modify the overall phenomenology of UHECRs significantly (e.g., Kachelrieß & Semikoz 2006; Blaksley & Parizot 2012). However, more subtle effects could occur and are addressed separately (see Sect. 5.4 for a partial discussion). Blaksley & Parizot (2012) have also shown that a distribution of maximum energies results in a modification of the effective source composition, but this effect is already covered by the range of compositions that we explore (see Sect. 3.2). In principle, one may also consider a cosmological evolution of the number density and/or power of the sources, e.g., following the star formation rate as a function of redshift or another law characterizing the UHECR sources. However, this is known to modify the constraints imposed by the data on the astrophysical models only (or mostly) by requiring a different intrinsic spectral index at the source (see, for instance, De Marco & Stanev 2005; Berezinsky et al. 2006; Allard et al. 2006), with no significant effect on the observed cosmic rays in the GZK energy range. This is due to the horizon effect, which results the UHECR flux in this energy range being completely dominated by nearby, hence almost contemporary sources (compared to the average source evolution timescale). Therefore, the anisotropy patterns should not be expected to depend significantly on the source evolution.

A more important effect on anisotropies should, however, be expected if one relaxed the assumption of continuous sources. A new set of models with impulsive sources (e.g., in scenarios where GRBs are the sources of UHECRs) could then be simulated. This is not explored in the present paper. We simply note that, leaving the other parameters unchanged, an impulsive source model should give rise to stronger anisotropies than the corresponding continuous source model does, because of the resulting greater effective source power (for the sources contributing to the observed flux at a given time) and the narrower range of energies observed at a given time from a given source (from energy-dependent diffusion effects), and thus the narrower range of angular deflections. Various aspects of transient sources have been discussed in the literature (e.g., Kalli et al. 2011; Murase & Takami 2009; Takami & Murase 2012).

3.2. Source composition and energy spectrum

The relative abundance of the various nuclei accelerated at the UHECR sources – simply referred to as the source composition – has a strong influence on the level of anisotropies that one can expect to measure. Obviously, scenarios in which UHECRs are dominated by protons are much more favorable than scenarios in which Fe nuclei are dominant: ceteris paribus, the smaller deflections of protons result in much tighter multiplets observable in the sky map, and facilitate the study of individual sources and their identification. In the absence of prior knowledge about the source composition, we explore a range of possibilities, choosing models with the main requirement that they reproduce the measured cosmic ray energy spectrum above the ankle, based on either the Pierre Auger Observatory or the High Resolution Fly’s Eye (HiRes) and Telescope Array (TA) data, which at first order appear to differ only by an overall shift in the energy scale (Dawson et al. 2013; Fukushima 2013).

The energy spectrum at the source is assumed to be a power law with a logarithmic index x, which is the same for all nuclear species. For each nucleus, i, we set a maximum energy, E_max,i, proportional to its charge, Z_i, so that E_max,i = Z_i × E_max,p, where E_max,p is the maximum proton energy (assumed to be the same in all sources). An exponential cut-off is then applied at E_max,i.

Here, we consider five different scenarios, referred to as

• 1.

  MC-high (mixed composition, with highE_max,p);

• 2.

  MC-4EeV (mixed composition with E_max,p = 4 EeV);

• 3.

  MC-15EeV (mixed composition with E_maxp = 15 EeV);

• 4.

  pure-p (pure proton model);

• 5.

  pure-Fe (pure iron model).

The “MC-high model” is the so-called mixed-composition model introduced in Allard et al. (2005), in which the UHECR source composition, below E_max,p, is assumed to be similar to that of low-energy Galactic cosmic rays. A good fit of the UHECR spectrum data is obtained by assuming a spectral index x = 2.3 and a maximum energy E_max = Z × 10^20.5 eV. It follows from this high value of the maximum proton energy that the composition is dominated by protons at all energies. Although the propagated spectrum expected for this model (Allard et al. 2005, 2007b,a) is compatible with the observations, the expected composition above 10¹⁹ eV does not appear to reproduce the Pierre Auger Observatory data concerning the average penetration depth of the induced atmospheric shower, ⟨ X_max ⟩, or its shower-to-shower fluctuation, rms(X_max), as a function of energy (Pierre Auger Collaboration 2010a). In principle, this model can, however, be reconciled with the data if the results concerning these observables are regarded as showing evidence of a change in the underlying hadronic physics rather than in UHECR composition.

The next two source models belong to the category of the so-called “low E_max models” described in Allard et al. (2008), which refers to mixed-composition models in which the protons do not reach the highest energies. In these models, protons are accelerated by the sources only up to a maximum energy that is lower than the GZK cut-off energy range. In such scenarios, the source composition above 10¹⁹ eV is gradually becoming heavier, with a dominant contribution of Fe nuclei above, say, 30 EeV. This appears to conform to the evolution of the UHECR composition suggested by the Pierre Auger Observatory data.

In the case of the mixed-composition “MC-4EeV model”, the adopted parameters are identical to those used in Allard (2012), where E_max,i = Z_i × 4 EeV, and the spectral index required to reproduce the observed energy spectra is relatively hard, namely x = 1.4. The abundance of heavy nuclei at the sources must be larger than in the MC-high model, to avoid an observable drop in the overall spectrum above E_max,p. However, the relative abundances of the heavy nuclei (Z ≥ 2) at the source are assumed to be the same. This model was shown to reproduce the composition trend drawn from the X_max observations made by the Pierre Auger Observatory relatively well.

An intermediate case is also considered, namely the “MC-15 EeV model”, in which the maximum energy of the protons is set equal to 15 EeV, so that E_max,i = Z_i × 15 EeV. While in the MC-4 EeV model the maximum proton energy is such that the abundances of C, N, O, and the other intermediate nuclei is very low above 50 EeV, these elements are still present with a significant abundance above 50 EeV in the source composition of the MC-15 EeV model.

Finally, we explore two “extreme”, but instructive models. The “pure-p model” corresponds to a scenario in which only protons are assumed to be accelerated at the source. This is considered for reference as an extreme case of light UHECR composition. The corresponding spectral index is x = 2.6 (see Berezinsky et al. 2006), and the maximum energy is taken as E_max = 10^20.5 EeV (or larger). From the point of view of the propagated composition at the highest energies, say above 50 EeV, the pure-p model is very similar to the MC-high model, which is dominated by protons at all observable energies. The general anisotropy features obtained in both cases are thus also very similar. Therefore, we do not show the results corresponding to the pure-p model in this paper separately.

In the “pure-Fe model”, only Fe nuclei are assumed to be accelerated at the source. Although this is not a realistic scenario from the astrophysical point of view, it is used for reference as an extreme case of heavy-composition models. The observed UHECR spectrum can be reproduced well within the pure-Fe scenario (e.g., Allard et al. 2008) using a source spectral index x = 2.3 and a maximum energy for Fe nuclei of E_max = 26 × 10^20.5 eV (or above), so that secondary protons, produced by photodissociation during propagation, can have energies up to ~10^20.5 eV. This model is characterized by a composition at the highest energies that is dominated by heavy nuclei (A > 40), but with the presence of a significant fraction (25−30%) of secondary protons, which are much less deflected by Galactic magnetic fields. This feature has interesting consequences for the anisotropy patterns, similar to what would be obtained from a composite astrophysical scenario in which a few sources provide a subdominant component of protons at high energy, in addition to a dominant heavy component. Therefore, this scenario is also explored to provide some hint of what a more complex, but probably more realistic scenario would imply (see the discussion in Sect. 5.4).

3.3. Source distribution and density

The anisotropy patterns in the UHECR sky maps also depend on the distribution of the sources in the nearby universe and on their spatial density.

In the absence of any clear indication about the nature of the sources, the most natural choice is to assume that they are distributed in a similar way as ordinary matter. The distribution of matter is known to be nonuniform in the nearby universe. In our simulations, the UHECR source distribution in tridimensional space (direction and distance) is derived from the distribution of galaxies in the 2MASS Redshift Survey catalog (2MRS, Huchra et al. 2012). More specifically, we use the initial survey with cuts in the near infrared magnitude mag and in Galactic latitudes | b | ≥ 10°. This catalog, which contains 20 860 galaxies, is linked to the Extragalactic Distance Database (EDD, Tully et al. 2009) to obtain the distance of nearby galaxies (~10–20 Mpc/h), for which the peculiar velocities dominate over the cosmic expansion, while the Hubble law in the linear regime is used to estimate the distances of more distant galaxies. Each galaxy in the resulting catalog is represented by a black dot in the distance-luminosity plane in Fig. 1 (left, top panel).

To compensate for the missing sources in the Galactic plane, we follow the filling method described in Crook et al. (2007), which consists in randomly populating the original sample to enhance the galaxy number density behind the Galactic plane, in a way that reflects the density observed just below and above it. As for the intrinsic luminosity of the additional galaxies, it is drawn according to the luminosity distribution function of the 2MRS sample. This luminosity function is derived from the catalog itself and is shown in Fig. 1 on the right. This procedure generates an additional set of 2094 galaxies and ensures the continuity of the structures across the plane. These additional galaxies are shown in the lower left hand panel in Fig. 1. The complete catalog is displayed in Fig. 2.

A magnitude-limited survey is affected by radial-selection effects, since at each distance D, only galaxies brighter than an intrinsic magnitude M₀(D) can be detected with apparent luminosity K_s. We use the curve M₀(D) corresponding to mag to build a volume-limited sample from a magnitude limited survey. To illustrate, all galaxies located inside the lower left hand rectangle in Fig. 1, with absolute magnitude greater than M_K = −23.75 and located within D = 100 Mpc, make a complete sample of galaxies.

Similarly, volume-limited samples can be obtained by selecting the galaxies in the lower-left corner of rectangles built in the same way on Fig. 1 (left) for different distances, with a corresponding cut at higher or lower luminosity. The resulting numbers of galaxies are shown in the middle plot of Fig. 1, as a function of maximum distance. As can be seen, the maximum number of galaxies in the various volume-limited samples is obtained for a cut at 100 Mpc, which by coincidence turns out to correspond to the so-called GZK sphere, related to the indicative horizon of 100 Mpc corresponding to ~90 EeV protons or Fe nuclei. It contains 6720 galaxies, 1267 of which have been added by the previously explained Galactic plane completion procedure.

In our simulations, we used this particular sample of galaxies as the seed catalog within which the UHECR sources are randomly chosen, with different values of the source density, n_s. As already indicated, the complete set of galaxies contains 6720 galaxies, which corresponds to a source density n_s = 1.6 × 10^-3 Mpc^-3. This catalog is then completed beyond 100 Mpc with galaxies isotropically distributed on the celestial sphere and whose magnitude is drawn randomly among the distribution of magnitudes of the seed catalog. To explore a UHECR scenario with a source density of 10^-5 Mpc^-3 (respectively 10^-4 Mpc^-3), we thus simply select randomly 1 galaxy out of 160 (respectively 1 out of 16) in the catalog.

We do not need to assume that the actual UHECR sources are necessarily among the 2MRS galaxies. All the anisotropy analyses that we perform are investigations of the intrinsic anisotropy of the simulated data, so the actual position of the sources in the sky is not relevant. Only the relative positions and the global angular/distance distribution is important. We thus simply assume that the overall distribution of the sources is similar to that of the galaxies. Moreover, for the density under study, the number of sources is low, and we do not expect the anisotropy signal to be driven by the spatial distribution of sources in 2MRS.

4. UHECR propagation and sky maps

4.1. General procedure

Our main goal has been to simulate realistic UHECR sky maps and to quantify their intrinsic anisotropies. For this, we computed the propagation of the UHECRs from their sources to the Earth and determine their energy, nuclear type, and arrival direction taking the relevant energy loss processes, the possible change of nuclear species, and the deflections in the extragalactic and Galactic magnetic fields into account.

We used the Monte Carlo code presented in Allard et al. (2005) to compute the energy losses and photodissociation processes in the extragalactic photon backgrounds (see Allard et al. 2008; Decerprit & Allard 2011, for a more detailed description). We also computed the 3D geometrical trajectories as influenced by the magnetic fields using the fast integration method described in detail in Globus et al. (2008), where a comparison with a full numerical integration is given. This allowed us to keep track of the time dependence (i.e., redshift dependence) of the energy losses without having to assume rectilinear transport.

The propagation is treated numerically in two separate steps. In the first step, we propagate the UHECR protons and nuclei in the extragalactic medium, following them in energy, mass, and geometrical spaces from their emission at a given source located at a distance D (and Galactic coordinates l and b), at a redshift/time z. This provides us with the propagated flux of the UHECRs injected by the whole set of sources as they enter our Galaxy, characterized by their energy and mass distribution, as well as by the apparent arrival directions.

The second step takes the deflections by the Galactic magnetic fields into account and relates the UHECR arrival directions on a fictitious sphere representing the boundary of the Galaxy to the observed arrival directions on Earth. This is done by inverting the relation between the different directions in the sky and the directions of the particles at the entrance of the Galaxy, as derived from the back propagation of negatively charged nuclei in the Galactic magnetic field, as explained below. The resulting “transfer function” of the Galaxy can then be applied to the extragalactic UHECR sky map to produce the desired sky map on Earth.

Finally, we analyse the anisotropy of the simulated data set by searching for significant excesses in the angular two-point correlation function. In the next sections, we describe the ingredients of the procedure in more detail and give the results in the next section.

4.2. Propagation in the extragalactic magnetic field

The extragalactic magnetic field (EGMF) is poorly known, and its spatial distribution, intensity, coherence length, time evolution, and origin are all uncertain. Observations imply the presence of μG fields in the core of large galaxy clusters. However, the spatial extension of these large field regions and their volume filling factor in the universe are difficult to evaluate. Efforts have been made to model local magnetic fields using simulations of large-scale structure formation that include an MHD treatment of the magnetic field evolution (see the pioneering studies by Dolag et al. 2002; Sigl et al. 2004; or more recent calculations by Das et al. 2008; Ryu et al. 2008, 2010; Donnert et al. 2009). Some of these simulations are constrained by the local density/velocity field to provide more realistic field configurations in the local universe. These simulations rely on different assumptions regarding the origin of the fields and the mechanisms involved in their growth. They are ultimately normalized to the values observed at the present epoch in the central regions of galaxy clusters (see the discussion in Kotera & Lemoine 2008). The outcome of the different simulations strongly differ. In particular, the volume filling factors for strong fields (above 1 nG, say) vary by several orders of magnitude from one simulation to the other. In contrast, an interesting simple alternative to complex hydrodynamical simulations, offering more freedom to test different models of magnetic field evolution with local density, has been proposed by Kotera & Lemoine (2008).

In view of the above-mentioned uncertainties, we used a simplified approach, assuming that the universe is filled with a purely turbulent, homogeneous magnetic field. Admittedly, such a topology of the EGMF is not realistic, but since our main purpose was to study the effect on the UHECR sky maps of the magnetic deflection in the extragalactic space, we simply investigated the impact of a magnetic blurring upstream of the Galaxy for different typical values of these deflections. The smallest impact corresponds to no magnetic field, while the largest impact would be obtained for high EGMF intensities of a few nG. Large localized magnetic fields, on the other hand, could be important if the volume filling factor is not too small, but this would mostly result in the apparition of fake secondary sources at the location of the magnetic cores, and thus produce a similar phenomenology, with only a higher apparent source density (Kotera & Lemoine 2008). Therefore, we simply simulated a uniform turbulent EGMF, with a method following the one described in Giacalone & Jokipii (1999). We assumed a Kolmogorov-like turbulence with a maximum scale λ_max = 1 Mpc and different magnetic field variance, $\sqrt{\mathrm{⟨}{\mathrm{B}}^{\mathrm{2}}\mathrm{⟩}}$, ranging from 0.1 to 3 nG.

4.3. The Galactic magnetic field model

To simulate the transport of the UHECRs in the Galaxy, we implemented a representative model of the Galactic magnetic field (GMF). We followed the modeling of Jansson & Farrar (2012a,b) who consider three types of magnetic structures: i) a coherent field with spatial scales of a few kpc; ii) an isotropic turbulent field with spatial scales of tens of pc; and iii) a striated field, which refers to an anisotropic turbulent field whose orientation is aligned with the large scale coherent field, but whose strength and sign vary on a small scale.

The large scale coherent field is modeled as the superposition of three separate components: a disk component, a halo component and an out-of-plane halo component. The disk component originates in the Brown et al. (2007) model where the magnetic field is concentrated in the plane and closely follows the spiral arms of the Galaxy. Several large scale reversals of the magnetic field occur along the Galactic radius. The disk field is symmetric with respect to the Galactic plane and transitions smoothly to a strictly azimuthal (toroidal) halo field at low vertical extent. This halo field decreases exponentially with scale height and takes different amplitudes below and above the plane. Finally, the out-of-plane halo component is inclined with respect to the Galactic plane, with a constant inclination far from the Galactic center and an almost perpendicular orientation closer to the Galactic axis. This so-called X-field is motivated by observations of X-shaped field structures in external galaxies (Krause et al. 2006; Krause 2007; Beck 2009). The large scale regular field model of Jansson & Farrar (2012a,b) – i.e., the sum of the disk field, the toroidal halo field and the X-field – is illustrated in Fig. 3.

For the random component of the GMF, we follow the modeling procedure of Giacalone & Jokipii (1999) and assume a Kolmogorov-like turbulence with a coherent length λ_c = 100 pc. Its intensity is given by the field variance, which we assume follows the magnitude of the regular component of the GMF, with an overall enhancement factor of 3. In other words, the turbulent field is typically rms magnitude 6 μG where the regular field is 2 μG. This scaling is assumed to be isotropic, so that magnetic turbulence associated with this random component is isotropic and spatially homogeneous on small scales.

The last component of the GMF model is the striated field, which is included after the recent works of Jaffe et al. (2010, 2011). Here, the field is either parallel or antiparallel to the regular field with a coherence length of 100 pc, and its magnitude follows the magnitude of the regular component, according to Jansson & Farrar (2012a,b).

4.4. Particle deflections in the Galactic magnetic field

To build the UHECR sky maps, we need to connect the arrival direction of the particles into the Galaxy, as resulting from the extragalactic propagation, with the direction in which they are eventually observed on Earth. This is done in a statistical way by applying the following procedure.

4.4.1. Trajectories and global deflections

First, we back-propagate a very large number of protons away from the Earth, until they reach a sphere of radius 50 kpc centered on the Galactic center, loosely considered as the “boundary” of the Galaxy, beyond which the influence of the GMF is negligible. More specifically, we propagate antiprotons with fixed energies between 10^17.5 eV and 10^20.5 eV, by steps of Δlog (E/ [ eV ] ) = 0.1, starting from the Earth in different directions. For each energies, the starting directions are regularly distributed over the celestial sphere using an HEALPix grid (Górski et al. 2005) with resolution parameter N_side = 1024, which corresponds to 12 582 912 directions, or a pixel size of ~3.5 arcmin. The spatial transport of the particles is then computed by simply integrating the equation of the trajectory governed by the Lorentz force (see Stanev 1997; Harari et al. 1999; Takami & Sato 2008; Giacinti et al. 2010b, for a discussion on cosmic ray propagation in the GMF).

In Fig. 4, we show a set of ten trajectories of back-propagated 5 EeV (anti)protons bent by the GMF. The distance traveled by the particles (along their trajectory, but measured in the Galactic frame) is always greater than the rectilinear distance from Earth, because of the deflections. This is shown in the lower right hand panel of Fig. 4. However, even when the deflections are relatively large, the difference does not exceed about 10%. Whatever the starting direction, the (anti)protons are found not to be confined by the GMF, and the residence time inside the Galaxy remains negligible compared to the energy loss time, so we can neglect their interactions with the local photon fields. This justifies that we only consider the Lorentz force when computing the particle transport in the Galaxy, and as a consequence, all UHECRs with the same magnetic rigidity behave in the same way. The trajectories of 5 EeV protons are thus also those of 30 EeV carbon nuclei, 40 EeV oxygen nuclei, 130 EeV iron nuclei, or any nucleus with a 5 EV rigidity.

The above back-propagation gives us, for each rigidity, a one-to-one relation between the ~12.6 million starting directions on the celestial sphere and the direction in which the corresponding (anti) UHECR leaves the Galaxy. In Fig. 5, we show the histogram of the resulting angular deflections for all these UHECRs for four different rigidities: 10^0.7 ≃ 5 EV (relevant for ≃130 EeV Fe nuclei), 10^1.2 ≃ 16 EV (relevant for ≃95 EeV C nuclei), 10^1.8 ≃ 63 EV, and 10^2.1 ≃ 130 EV.

As expected, the deflections are much larger at low rigidities, and UHECRs with intermediate rigidities experience a wide range of deflections, depending on the arrival directions. The difference between the top left hand panel and the bottom right hand panel illustrates the difference between a proton primary and an iron nucleus primary at the highest energies. Obviously, direct pointing astronomy seems inaccessible if the UHECRs are dominated by Fe nuclei. However, we recall here that direct source pointing should not be the only goal of UHECR astronomy, and the study of anisotropy patterns can provide important information about the UHECR sources. In addition, a small number of protons (or low rigidity nuclei) may lead to the appearance of small angular scale multiplets if the statistics are large enough to allow detecting a few of them. Moreover, even in the most unfavorable case where all UHECRs above ~80 EeV are Fe nuclei, a noticeable number of them are found to experience deflections smaller than the typical angular separation between sources. To illustrate this point in a more quantitative way, we show the fraction of arrival directions corresponding to UHECR deflections smaller than 10, 20, 30, 40, 50, and 60 degrees in Fig. 6, as a function of rigidity (translated into Fe nuclei energy on the top x axis). Finally, we note that some knowledge of the regular component of the GMF can in principle be used to correct for the non-random part of the UHECR deflections.

For comparison in Fig. 7, we show the deflection fractions of protons with energies between 1 and 300 EeV. At 50 EeV, 20% of the protons are deflected by more than ten degrees. This fraction drops to less than 3% above 100 EeV.

While the distributions illustrated by Figs. 5–7 mix all the arrival directions together, more information about the UHECR deflection patterns across the GMF can be obtained by looking separately at different pixels on the sky map.

4.4.2. Backward and forward deflection maps

Starting from the above one-to-one relation between the starting direction of the back-propagated particles and their direction out of the Galaxy, we then use a coarser sampling of the celestial sphere, choosing a HEALPix resolution parameter N_side = 64. This defines 49 152 pixels of slightly less than 1deg² evenly distributed over the sky, each of which contains 256 of the original directions on the fainter grid. Thus each direction on the sky (with a resolution of ~1deg) is now linked with 256 directions at the boundary of the Galaxy, which are in effect the arrival directions of UHECRs with the rigidity under consideration that would be observed on Earth in that direction (with the assumed GMF).

For each pixel in the coarse grid (observed direction), we computed the average angular deflection, i.e. the mean of the 256 angles between the incoming directions at the entrance of the Galaxy and the observed direction. The result is shown in Fig. 8, where we plot the map of these mean deflections in a color code for the same four rigidities as above. Although the patterns show similar shapes, which are associated with the structure of the regular component of the GMF, the color code spans different ranges in each map to follow the global reduction of the deflection with increasing rigidity.

It is interesting to note that the range of deflections is generally quite broad over the celestial sphere. UHECRs observed in some directions are on average much less deflected than in others. Particles observed in a large circle of ~15° radius around the Galactic center are notably much more likely to have been deflected by a large amount than particles observed toward anti-center longitudes, especially in the northern (Galactic) hemisphere. This strong contrast between observing directions is partly reflected in the distribution of the angular deflections of Fig. 5, where the top panels show two wide, but distinct peaks. But most importantly, it can in principle be exploited to perform refined anisotropy analyses attributing different weights to different regions, based on some prior knowledge of the relative deflection amplitude. This is not attempted here.

The above deflection maps may be misleading, however, since they give information about the average deflection of the UHECRs observed in different directions, but not about the UHECRs coming from sources located in these directions. For the same reason, the results of the back-propagation of charged particle cannot be exploited to build simulated sky maps until an inversion is done to relate the arrival directions of cosmic rays at the entrance of the Galaxy (from their particular extragalactic sources) to the actual directions in which they are observed on Earth. We may call a “Galactic pixel” a pixel in the sky map where a back-propagated particle exit the Galaxy, or where a forward-propagated particle enters the Galaxy. Likewise, an “Earth pixel” is a pixel in the sky map where a back-propagated particle starts its trajectory or where a forward-propagated particle is eventually observed.

The inversion is done by simply keeping track, for each UHECR incoming direction (i.e., for each Galactic pixel), of the different Earth pixels in the direction where back-propagated particles were initially sent to exit the Galaxy in that pixel. The number of Earth pixels related to a given Galactic pixel is not known a priori. On average, 256 pixels on the fine grid are associated with the Galactic pixels on the coarse grid. This is also, of course, the number of Earth pixels that would be associated with each Galactic pixel if there were no deflection at all. However, some directions turn out to be more likely fed by back-propagated particles than others, because the GMF can focus back-propagated particles from a wide range of directions into a smaller solid angle, or conversely. This is illustrated in Fig. 9, which shows the magnification factors for each Galactic pixel, defined as the number of Earth pixels associated with that pixel divided by 256. For instance, a magnification factor of 2 indicates that a source in that direction will contribute twice as much flux to the UHECRs observed on Earth as if there were no deflection (or as the average of what could be expected if it were anywhere in the sky).

There is a different magnification map for each particle rigidity and, as can be seen, the level of magnification can vary by large amounts even between two relatively nearby source directions. This is the well known phenomenon of the so-called caustic curves (see Harari et al. 1999, for instance), which are singular mathematical lines where the wave front of a radial stream of particles bent by the (regular-only) magnetic field would be tangent to the line of sight. Strongly contrasted caustics are found to appear only for intermediate rigidities in the energy range of interest. Indeed, for low rigidities, the particles are strongly deflected along any direction and, in the limit of very large deflections, an essentially isotropic flux is produced and the magnification tends to one in all directions. Conversely, at very high rigidity, the deflections become very small in any direction, and the particle transport across the GMF tends to the trivial one-to-one relation between Earth pixels and identical Galactic pixels. At intermediate rigidities, complex structures can be observed, with regions of large magnification next to regions of large demagnification.

Quantitatively, one can see from the top panel of Fig. 9 that, for some source locations, the magnification factor can reach values as high as 10 to 18 for Fe nuclei around 130 EeV, while we appear to be essentially blind to other regions of the sky. The same is still true at a rigidity of 16 EV, with magnification factors up to 9.6. The magnification factors do not exceed the value of three or four at rigidities over 60 EV, but regions of the sky, notably just below the Galactic plane, appear to be strongly demagnified, down to almost complete invisibility, even for protons that are very little deflected in this energy range.

It is worth noting that there is no contradiction with the so-called Liouville theorem, which indicates that an incoming isotropic flux must remain isotropic, regardless of the intervening magnetic field. This simply means that dead zones (or angles) from which cosmic rays are deflected away and which never reach the Earth are compensated exactly by focused zones from where incoming cosmic rays are deflected into the apparent direction of the dead zone. Here, the sources are discrete, and some positions can be partially or completely hidden to us by the GMF, while others can be magnified to higher apparent luminosities. In principle, this can modify the source statistics by potentially large amounts, depending on the range of rigidities of the UHECR particles. In addition, since the magnification of a source in a given direction depends so much on the rigidity, the magnification factor can be very different for different types of particles at a given energy. In particular, a given source location can result in a strong magnification (or demagnification) of the Fe component in a specific energy range, while the proton component will be unaffected. This may result in noticeable modifications of the composition and spectrum of individual sources or of separate regions of the sky. These interesting effects and the associated constraints that may be derived from them will be studied in a statistical way in a forthcoming paper.

In addition to the magnification factors, the inverted (i.e. forward) relation between Galactic pixels and Earth pixels allows us to determine where the UHECRs entering the Galaxy in a given direction will be observed on Earth. From the set of Earth pixels associated with a given Galactic pixel, it is straightforward to compute the average deflection experienced by the UHECRs coming from a source in that direction, as a function of rigidity. This is shown synthetically in Fig. 10 in the form of so-called “forward deflection maps”, where we plot the mean of the set of deflections of UHECRs not observed in a given direction (as in the “backward deflection maps” of Fig. 8), but initially coming from that direction. This mean is calculated from the set of Earth pixels associated with that direction, which contains on average 256 pixels but can be much more or much less numerous depending on the magnification factor. Some particularly blind source directions turned out to be associated with no Earth pixels at all, making it impossible to derive a mean deflection. The corresponding pixels have been represented in gray in Fig. 10.

As mentioned above for the backward deflection maps, an important feature of the forward deflection maps is the strong contrast between the observational situation for different source directions. Even at the lowest rigidity represented here (top panel), where deflections are on average very large, some significant parts of the sky give rise to much smaller deflections, represented in blue in the plots. The prior knowledge of these regions (which derives directly from the knowledge of the GMF) should help in deriving meaningful constraints from the distribution of UHECR arrival directions and anisotropy patterns.

Finally, it is interesting to note that large deflections usually correspond to either randomization, in which case the magnification tends to be close to unity, or to obscuration, in which case the magnification factor can become much lower than 1, or even tend to zero (e.g., toward the Galactic center, behind which a source has very low probably to be visible at UHE). On the other hand, large magnification usually require an ordered variation in the deflections occurring over a range of nearby directions, which can coherently extend the solid angle “feeding” a given direction.

Of course, the exact patterns observed on these various maps (backward, forward, and magnification maps) depend on the GMF model used in the propagation code, which is unlikely to be correct all across the Galaxy. However, we may hope that the assumed GMF model is sufficiently representative of the actual GMF for the above results to give a reasonable idea both of the typical average deflections and standard deviation values and of their range of variations over the sky.

4.5. Generation of UHECR sky maps

4.5.1. General procedure

The final step in the propagation procedure is to build the simulated sky maps, each of which represents a particular set of UHECR events distributed over the sky, as could be observed by a given experiment. For this, we simply put together the above elements.

First, we select an astrophysical scenario; i.e., we choose one of the five generic composition/spectrum models described in Sect. 3.2 and summarized in Table 1, and assume a given source density, n_s. We then generate a particular realization of this scenario by selecting the location of the sources through a random draw in the source catalog with density n_s, as explained in Sect. 3.3. A total of 500 different realizations are simulated for each astrophysical scenario, in order to explore the so-called cosmic variance of the models, i.e., the range of sky map properties that can be expected within a given scenario, depending on the contingent distribution of the actual sources currently active in the local universe.

For each realization, we build different sky maps, depending on the intended observatory (determined by its coverage map), and on event statistics (determined by the total exposure of the experiment). In the current paper, we chose either the partial sky coverage of the Pierre Auger Observatory or the almost uniform sky coverage of JEM-EUSO with an exposure of 300 000 km² sr yr, as discussed in Sect. 1.

The UHECR particles are generated one by one, with their own energy and nuclear type, according to the source spectrum and source composition of the scenario under investigation. Each particle is propagated in the EGMF with the Monte Carlo code described above. We then applied the magnification factor appropriate to the resulting direction of the UHECR as it enters the Galaxy, as derived in Sect. 4.4.2. For this, we normalize the magnification map to the maximum magnification factor at all energies and apply a standard acceptance/rejection method. To determine the observed direction of the UHECRs on Earth, we randomly chose between the various Earth pixels associated with the incoming direction (i.e., Galactic pixel). The next step consists in applying the coverage map of the experiment, i.e., accepting/rejecting the events according to the normalized exposure in the relevant arrival direction. For the JEM-EUSO-like detector, we applied an additional acceptance/rejection procedure to account for the detection efficiency as a function of energy. For this, we used the efficiency curve computed for JEM-EUSO, as given in Adams et al. (2013). Finally, we applied an error on the energy and direction to reflect the experimental uncertainty on the reconstructed shower parameters. For JEM-EUSO, we used a simplified and conservative Gaussian energy resolution uncertainty of 30% for all the UHECR events. The angular resolution is also assumed to be Gaussian with a width of two degrees, but given the patterns and amplitudes of the deflections for the different models, changing this parameter did not appear to have any significant impact on the results.

The above procedure was applied to each UHECR, one after the other, until the intended statistics are collected for the detector under consideration. The resulting sky map is the final output of the simulation, from which a systematic search for anisotropy can be performed, as discussed in Sect. 5.1.

4.5.2. Sky map statistics

It is important to note that the total number of events to be observed by a given observatory is not known a priori. The main source of uncertainty on the UHECR flux resides in the so-called energy scale of the measured UHECR spectrum. Both the Pierre Auger Observatory energy spectrum and the TA energy spectrum suffer from systematic uncertainties on the reconstructed energy of the UHECR events. While a joint working group has proposed to build a fiducial energy spectrum by rescaling the Pierre Auger Observatory energy scale upward and the TA energy scale downward (Dawson et al. 2013; Matthews 2013), the exact flux remains uncertain. In this paper, we considered the two assumptions on the energy scale separately (i.e., we did not apply any rescaling), which result in two different assumptions for the UHECR flux as a function of energy, hence two different statistics expected above some energy threshold, for a given choice of the total exposure.

Another source of uncertainty in the expected number of events is the absence of clear knowledge of the shape of the UHECR spectrum, which requires larger statistics to be known precisely and may be different in different regions of the sky. As a matter of fact, detecting anisotropies in the UHECR arrival directions above a given energy is equivalent to detecting a different energy spectrum in different directions. The results discussed below show that significant anisotropies should be expected at the highest energies, regardless of the assumed astrophysical scenario, so we cannot use the current knowledge of the spectrum in a limited region of the sky (even barring its imperfection) to predict the number of events that an all-sky coverage experiment should detect. The uncertainty associated with the current poor knowledge of the shape of the spectrum is, however, much smaller than what is associated with the energy scale. For our present purpose, we simply assume that either the Pierre Auger Observatory flux or the TA flux hold over the whole sky, and derive a fiducial spectrum by averaging the spectra obtained by fitting the currently available data with our different models. This fiducial spectrum is then used to determine the expected numbers of events for a given detector. Two reference spectra are thus built, one for each choice of energy scale. From these, we determined the following statistics for the JEM-EUSO-like detector with the quoted total exposure. In the case of the Pierre Auger Observatory energy scale assumption, we expect 1100, 250, and 100 events, respectively, above 50 EeV, 80 EeV and 100 EeV (implementing the energy detection efficiency of JEM-EUSO (Adams et al. 2013). In the case of the TA energy scale assumption, we expect 2100, 580, and 260 events above the same energies. Although model-dependent, these numbers represent the best-guess limits on the number of events that can be extrapolated from the current knowledge of the spectra measured with low statistics, partial sky-coverage ground observatories, within the framework of the astrophysical scenarios investigated here.

Finally, in addition to these fiducial statistics for a future observatory similar to the JEM-EUSO project, we also consider the current Pierre Auger Observatory statistics as a reference point for selecting astrophysical models that appear compatible with the current data, as far as anisotropy is concerned (see below). For this, we apply the Pierre Auger Observatory coverage map and accumulate 84 UHECRs above 55 EeV, which corresponds to the statistics gathered by the Pierre Auger Observatory up to June 2011 (according to Kampert 2011, in which the last search for small scale anisotropy above 55 EeV with the Pierre Auger Observatory is reported).

4.5.3. Reading the sky maps

An example of a set of sky maps is shown in Fig. 11. This is the result of a typical simulation, corresponding to a particular realization of a mixed-composition model with a source density of n_s = 10^-5 Mpc^-3 and a maximum proton energy of E_max,p = 15 EeV (MC-15 EeV model). The map in the top panel is the Pierre Auger Observatory-like reference map, showing the arrival direction of 84 UHECR events above 55 EeV. The map is shown in Galactic coordinates, and the wide region without events in the left and upper right parts of the map are regions of the sky inaccessible to the detector.

The symbols and color codes obey the following rules:

• the shape of the symbols representing the events give an indication of the mass of the associated UHECR: polygons with larger numbers of sides correspond to heavier nuclei, as indicated on the map, and protons are shown as circles;

• the size of the symbols is proportional to the particle energy: larger symbols correspond to higher energy particles;

• events shown with the same color correspond to UHECRs coming from the same source; however, only the most intense sources (by decreasing multiplicity and provided that they contribute at least three events and 1% of the total flux) are shown with a separate color; all the other events are shown in black;

• the colored stars correspond to the real location of the sources of the events sharing the same color.

As can be seen in the top panel of Fig. 11, in this particular example only three sources have a multiplicity higher than four in the map corresponding to the Pierre Auger Observatory statistics of reference. The most intense one is shown in red, with eight events out of the 84 events recorded. The second source is in blue, with seven events, and the third source is in green, with four events in the field of view. It is interesting to note that the source of these 4 events color-coded in green is actually far away in the unobservable part of the sky. The distance to the color-highlighted sources is also indicated in the map. In this case, the “green source” is the closest, most luminous source in the sky, at 23 Mpc. The four events observed from this source in the southern hemisphere sky consist of one low-mass (square symbol), one intermediate mass (pentagon), and two heavy (hexagon) nuclei, deflected in the Pierre Auger Observatory field of view by the GMF. The rightmost event has the lowest energy, as indicated by its smaller size.

No obvious clustering of events is visible in the map, which is compatible with the absence of any clear small-scale anisotropy in the Pierre Auger Observatory data. In this model, this is mostly due to the low value of the maximum proton energy assumed, namely 15 EeV, which results in the dominant presence of heavy nuclei in the energy range under consideration, as can be checked directly on the map (polygonal symbols). Two protons (circles) can nevertheless be seen in red, very close to their actual source, represented by the red star. With these statistics, such a doublet cannot be used to pin point a source on the sky, since the arrival direction coincidence could be a simple chance coincidence. As a matter of fact, a few other doublets can be observed in the sky map, including an almost perfect (but random) coincidence of two heavy nuclei (hexagons) coming from different sources. Three other isolated protons can be seen in the map. They may be close to their actual source, but the small statistics do not allow to identify these sources. Larger statistics are needed to assess the presence of real multiplets with reliable statistical significance.

The three subsequent panels in the same Fig. 11 show the sky maps obtained from the same realization of this astrophysical model (same source composition/spectrum/density scenario and same sources), assuming a uniform full-sky coverage with the JEM-EUSO-like exposure and detection efficiency, with energy thresholds at 50 EeV, 80 EeV, and 100 EeV, respectively. The three maps illustrate the GZK horizon effect, by which the distant sources contributing to the observed flux become less and less numerous as the energy increases. Even though the nearby, high-luminosity sources are also present, and even dominant at lower energies, the contribution of a very large number of sources distributed more or less uniformly makes it difficult to isolate the brightest sources. Even though a standard test of anisotropy should easily reveal the presence of significant excess in the maps (see below), it may not be as easy to derive meaningful astrophysical information from the map drawn with an energy cut-off of 50 EeV, as from the map drawn with an energy cut-off of 100 EeV (bottom panel).

At the highest energies, the dominant source appears to account for 40 of the 100 events (in agreement with the general results about UHECR source statistics reported in Blaksley et al. 2013). The corresponding source (red-colored events) is so-to-say self-isolated on the sky, because the more distant sources are cut off by the GZK effect. One should not be misled by the color code, however. In practice, of course, we will have no way, a priori, to distinguish the events coming from a given source. The green and purple sources appear much more mixed together and can not be so easily isolated. This was to be expected anyway, since their angular separation is close to the typical deflection of the individual events (here dominated by heavy nuclei), and their distance is roughly similar (33–34 Mpc), so that their apparent luminosity is almost identical and the GZK effect operates in the same way for both of them. This is a standard case of source confusion. The possibility to isolate (to a large extent) the dominant source in the sky at 100 EeV is nevertheless a common feature of most of the models and realizations that we have generated. It is worth pointing out here that each sky map is built from a different run; i.e., we draw a whole new set of events for each map. Therefore, even though the source realization of the model under study is the same for all the maps of a given figure (which means that the sources of the catalog selected as actual UHECR sources are the same), the events seen above 100 EeV in the map built with an energy threshold of 80 EeV, say, are not the same as those shown in the map built with a threshold at 100 EeV, which is in all respects a different map, from a different data set. In other words, the data sets for the different maps are all totally independent. For this reason, a direct comparison of the maps at different thresholds cannot be done here visually, but can of course be carried out by applying an energy threshold on the lower energy maps.

While the sky maps are sometimes useful to guide the eye, the best way to obtain definite and objective information about the UHECR sky is usually to perform unambiguous anisotropy analysis. This is what we did to analyse the thousands of sky maps generated by the above procedure in a systematic way, as explained in the next section.

5. Results and discussion

Our goal is to determine whether a future experiment with a total exposure of 300 000 km² sr yr could observe significant anisotropies in the UHECR arrival directions. This depends not only on the general astrophysical model assumed, but also on its particular realization, i.e., on the specific distribution of the sources, which would happen, within this model, to contribute to the observed flux of cosmic rays in our specific location in the universe, during the time of observation. Indeed, even for a given assumption about the source composition, spectrum, and density, a relativity wide range of situations could be encountered, which exhibit either very strong, moderate, or low anisotropies depending on this particular source distribution.

In the present study, we explore the anisotropy properties associated with the different scenarios in a statistical manner. To this end, we simulate 500 different realizations of each scenario and determine the probability that a given realization, chosen randomly, gives rise to a significant anisotropy.

5.1. Intrinsic anisotropy searches

There are many ways to search for anisotropies in a data set (e.g., Abbasi et al. 2004, 2006; Pierre Auger Collaboration 2008a, 2012; Abu-Zayyad et al. 2012b, 2013b; Sommers 2001; Gorbunov & Troitsky 2005; Harari et al. 2006; Kashti & Waxman 2008; Ave et al. 2009; Golup et al. 2009; Berlind et al. 2010; Giacinti et al. 2010a; Oikonomou et al. 2013), and for any given sky map, various tests can be performed. This includes searches for specific correlations with known astrophysical sources or for any type of potentially meaningful pattern noticed a posteriori on a partial data set, which can then be tested on subsequent, independent data sets as the sky map is being built over time. It is thus not possible here to do an exhaustive search for anisotropies, and most of the tests based on correlations with other catalogs of sources would mostly be irrelevant in the case of our simulated sky maps, not mentioning that the knowledge of the magnetic field is still incomplete at the moment. We thus chose to limit our studies to intrinsic anisotropies, as can be revealed by analysing the angular autocorrelation of the UHECR arrival directions, through localized excesses in the number of events or through the commonly used two-point correlation function.

The results shown here are conservative in this respect, in the sense that they do not exploit the whole arsenal of analysis tools available, but consider only the most standard ones, which are largely independent of the assumptions made on the GMF (see the related comment in Sect. 4.4). In addition, these tests do not depend much on the actual choice of the sources, but rather on their overall properties and density in the nearby universe. Likewise, in our anisotropy searches, we did not make any use of the potential correlation between the arrival direction of the UHECRs and their energy, which should display, on average, some coherent patterns in the case of UHECRs originating in the same sources (Harari et al. 2006; Kachelrieß et al. 2007). Such correlations, which should open new possibilities for identifying sources or constraining their number and locations, will be studied separately.

Here, we present the results obtained with one of the most widely used statistical tests for measuring the departure of a given data set from isotropy, based on the so-called two-point correlation function. This function characterizes the autocorrelation of the UHECR arrival directions by simply giving the number of pairs of UHECR events separated by less than a given angle on the celestial sphere, as a function of that angle. The test consists in comparing the observed numbers of pairs of events on any angular scale with the numbers of pairs present on the same angular scale in datasets of the same size built randomly from an isotropic distribution.

To ensure good statistical power for this test, we computed 10⁶ different realizations of an isotropic sky for each number of events considered (corresponding to different exposures and energy thresholds), implementing the appropriate coverage map, either that of the Pierre Auger Observatory or the uniform full-sky coverage relevant to JEM-EUSO. Each of these realizations has its own two-point correlation function, and the whole set of realizations gives us not only the average number of pairs of events on each angular scale, but also the distribution of the numbers of pairs that can be expected for an isotropic sky. From this distribution, we can determine the probability that a given sky map could be built from an underlying isotropic distribution by simply counting the fraction of isotropic realizations that lie further away from the average distribution than the data set under study.

On the whole set of sky maps, we also performed a standard blind search test for localized excesses in some directions of the sky, as a function of both energy and angular window size, and analysed the corresponding significance of the detected anisotropy using the Li & Ma (1983) statistics. The results obtained were on average very similar to those derived from the two-point correlation function. Even though some particular sky maps appeared to show larger departures from isotropic expectations with one test rather than with the other, no additional information could be derived once averaged over the different realizations. Therefore, we only show here the results obtained with the autocorrelation function.

Obviously, scenarios in which Fe nuclei dominate are less likely to generate strong anisotropies, notably on small angular scales, than scenarios in which protons dominate. The deflection angle of individual particles is, however, not the only important parameter. The source density also has a strong influence on the possibility to detect significant anisotropies. In the case of a low source density, say n_s = 10^-6 Mpc^-3, very few sources contribute to the observed flux, and multiplets with high multiplicity are bound to be detected by observatories with increased exposure, making it much more likely for clusters of events that are incompatible with an isotropic flux to be observed (Takami et al. 2012). Conversely, a very high source density will result in much lower multiplicities, even for the most luminous sources, and in addition a smaller angular separation between sources on average, making source confusion much more likely. To be definite, in the following we loosely refer to the case where n_s = 10^-6 Mpc^-3 as the low-density case, to the case where n_s = 10^-5 Mpc^-3 as the intermediate-density case, and to the case where n_s = 10^-4 Mpc^-3 as the high-density case.

5.2. Mixed-composition models with low E_max

The initial hope behind the intense observational efforts developed in the past decades to detect UHECRs was to identify their sources by direct pointing, in an energy range where the deflections by intervening magnetic fields would be small enough to let us identify tight clusters of events associated with individual sources, right behind their arrival directions. Such an ideal situation has not occured so far. A reason for this may be that the deflections are larger than initially anticipated, possibly because a dominant fraction of the highest energy cosmic rays are not protons, but heavier nuclei, with a lower rigidity. If this is the case, then the absence of any clear anisotropy detected by the current detectors is quite easy to understand. The key question is now whether a new generation of detectors, increasing the exposure by a factor of ten or so, could change the situation enough. In this section, we study the case of mixed composition models, in which the maximum energy of the protons in the sources is lower than the GZK energy scale, so that the highest energy particles are dominated by heavy nuclei and, most particularly, by Fe nuclei.

5.2.1. High source density: n_s = 10^-4 Mpc^-3

In Fig. 12, we show the two-point correlation functions obtained with the mixed-composition, low E_max,p model, MC-15 EeV, assuming a source density of 10^-4 Mpc^-3. Curves with different colors correspond to different statistics, and for each angular scale, the error bars contain 90% of the 500 realizations of that particular astrophysical scenario. These curves are interesting only in comparison with the isotropic expectations, which are shown by the shaded areas of the same color. These areas contain 90% of the million isotropic realizations.

The two lowest curves, in blue, correspond to sky maps simulated with the Pierre Auger Observatory statistics. As can be seen, a large fraction of the realizations are compatible with isotropy, and the average two-point correlation function for this model lies only marginally away from the isotropic expectations, especially if the Pierre Auger Observatory energy scale (light blue) is assumed. This confirms that the MC-15EeV scenario is fully compatible with the current data, as far as the anisotropies based on the autocorrelation of UHECRs are concerned.

The other six curves, however, show that with the statistics of JEM-EUSO, this scenario would produce significant anisotropies for most realizations, if not all. This is true for the sky maps built with a cut-off at 50 EeV (red curves), 80 EeV (beige curves), or 100 EeV (green curves), regardless of the assumption on the energy scale (lighter tone for the Pierre Auger Observatory scale, darker tone for the TA scale): the lower limit of the error bars hardly touches the shaded areas corresponding to the isotropic expectations, for most of the angular scales.

Another way to look at these results is proposed in Fig. 13, where we show the number of realizations of the same scenario that are further away than 2σ, 3σ, or 4σ from the average isotropic expectation. Since the distributions are not necessarily Gaussian, what we mean by n-σ is actually not the number of standard deviations away from the average of the distribution, but the distance away from the average that corresponds, respectively, to 95.5%, 99.7%, and 99.994% of the isotropic realizations. In other words, the non-Gaussian tails are taken into account to provide real probabilities rather than standard deviations.

The various curves in Fig. 13 confirm, in a more quantitative way, that the scenario under investigation is generally compatible with the Pierre Auger Observatory constraints on anisotropy, since fewer than 20% of the realizations show an anisotropy with a significance of 3σ, and fewer than 40% show an anisotropy as weak as 2σ (for the two-point correlation function). Assuming the TA energy scale instead of the Pierre Auger Observatory one increases the number of realizations showing some anisotropy, because the considered events are more energetic and thus less deflected and because the corresponding horizon distance is somewhat reduced. However, more than 50% of the realizations are still found not to display any significant autocorrelation.

The situation with the JEM-EUSO statistics is very different, since more than 80% of the realizations display at least a 3σ anisotropy at all energies. If one assumes the TA energy scale (dark colors), one even finds that 90% of the realizations can be declared anisotropic with a 4σ significance (i.e., a 6.3 × 10^-5 chance probability). Even with the Pierre Auger Observatory energy scale, more than 50% of the realizations show 4σ-significance anisotropy at 100 EeV, despite the lower statistics. This fraction increases to ~70% at 80 EeV, and up to more than 80% at 50 EeV. By contrast, the Pierre Auger Observatory is expected to detect a 4σ anisotropy in fewer than 10% of the realizations of such a scenario.

As can be seen, the angular scale where most realizations display the strongest anisotropy in this astrophysical scenario is relatively large, around 15–25 degrees. This is not surprising, given the predominance of high-Z nuclei among the UHECRs.

In Fig. 14, we propose yet another way to look at the results by showing the number of realizations that display an anisotropy stronger than a given significance, as a function of that significance (translated into a number of sigmas, as explained above). Since the significance depends on the angular scale, we chose for this plot the angular scale giving the maximum significance. In principle, one should penalize the resulting probability for searching at that particular angular scale (i.e., we should marginalize on angular scales). However, as clearly shown by Fig. 13, the angular scale where the maximum significance occurs is essentially always the same for the model under study (which is also true for the other models), and the curves are very flat over a wide range of angles. The departure from anisotropy could thus be searched for a priori in this angular range by fixing the angular scale before any trial, say at 20°, in which case no penalty factor should be applied. In conclusion, for each given astrophysical scenario tested, the penalization to be applied to Fig. 14 and to similar plots for other scenarios should remain limited.

The curves in Fig. 14 confirm the previous conclusion: while the Pierre Auger Observatory statistics are too low to allow the detection of any significant anisotropy in this MC-15EeV scenario, except in some “lucky” realizations representing only a small fraction of the possible skies, increasing the exposure up to that of JEM-EUSO will almost certainly lead to the detection of significant anisotropies. This is particularly true if the TA energy scale holds, as shown by the three darker lines, almost superimposed close to the 100% probability line. The JEM-EUSO statistics thus appear appropriate to efficiently constraining such a scenario.

It is worth noting that this scenario is one of the worst possible scenarios as far as the detection of anisotropies is concerned, since it implies that the high-energy particles are mostly heavy nuclei (thus experiencing large deflections), and the source density is high. Nevertheless, a detector able to reach the statistics of JEM-EUSO would be sufficient to detect significant anisotropies in essentially all the realizations of this scenario.

We also studied the effect of the EGMF. As discussed in Sect. 4.2, the most extreme situation may be given by a uniformly distributed, random magnetic field with a root-mean-square amplitude of 3 nG. Figure 15 shows how the results of Fig. 14 are modified in this case. As expected, the number of realizations that display significant anisotropy is somewhat reduced, compared to the case with a lower magnetic field. This is due to the larger dispersion of the UHECR arrival directions at the entrance of the Galaxy, which in turn results in a more widely spread distribution of arrival directions on Earth. However, the conclusions essentially remain the same. More than 90% of the realizations of the scenario under study display anisotropy with a significance greater than 4σ in the case of the TA energy scale, and almost 80% of them do so above 50 EeV in the case of the Pierre Auger Observatory energy scale.

5.2.2. Low source density: n_s = 10^-6 Mpc^-3.

The same plots as presented in the previous section can be built for all the scenarios. In Fig. 16, we show the results corresponding to the same low E_max,p model, MC-15 EeV, but with a low source density: n_s = 10^-6 Mpc³. As expected, the different realizations display more significant anisotropies overall, as a result of the smaller number of sources contributing at the highest energies, and thus of the higher multiplicity of the most intense sources.

The Pierre Auger Observatory data are still compatible with such a scenario, since a majority of the realizations (~60%) are not expected to have a two-point correlation function that departs from the isotropic expectation by more than 3σ, and ~40% of the realizations remain within 2σ of the isotropic expectation. It is interesting to note, however, that confirmation of the TA energy scale would create some tension between this astrophysical scenario and the Pierre Auger Observatory data, since only a bit more than one third of the realizations would then be compatible with a two-point correlation function less than 3σ away from the isotropic expectation.

Turning to the JEM-EUSO-like exposure, it clearly appears that such a scenario could be easily and severely constrained with the resulting statistics, if no significant anisotropy were detected. Essentially 100% of the realizations display anisotropies with a significance over 4σ, at all energies considered.

5.2.3. Comparison between the MC-15 EeV and MC-4 EeV models

The results shown above correspond to the MC-15 EeV model, i.e. to the intermediate value of the maximum proton energy achieved in the source, E_max,p = 15 EeV. Very similar conclusions can be reached for the lower value, E_max,p = 4 EeV, i.e. for the MC-4 EeV model. For comparison, we show in Figs. 17 and 18 the number of realizations of the MC-4 EeV model as a function of the significance of the measured anisotropy for the two extreme source densities, low and high.

For the lowest density, again, all realizations display a very significant anisotropy, regardless of the energy and angular scale. At the highest source density, the significance is reduced, as expected, even though a large majority of the realizations would be very significantly anisotropic, especially in the case of the TA energy scale.

However, an interesting feature distinguishes both scenarios. The main difference between the MC-4 EeV and MC-15 EeV models is the presence of relatively light nuclei up to higher energies in the latter case, due to the higher energy of the proton cut-off. In particular, C, N, and O nuclei, which are rather abundant in the interstellar medium, and thus in the assumed cosmic ray source composition, are accelerated in the sources up to respectively six, seven and eight times E_max,p. These intermediate nuclei are thus present up to ~100 EeV in the MC-15 EeV scenarios, while they disappear at ~30 EeV in the MC-4 EeV scenarios.

As an example, we show typical sky maps in Figs. 19 and 20 built from the MC-4 EeV and MC-15 EeV models at the intermediate source density, assuming the JEM-EUSO statistics with an energy threshold at 80 EeV (left panels) and 100 EeV (right panels).

The particular realizations chosen for these sky maps are right in the middle of the distribution of the 500 realizations of the corresponding scenario relatively to the significance of their anisotropy (as measured with the two-point correlation function). In other words, half of the realizations show stronger – i.e. more significant – anisotropies, while the other half show weaker anisotropies. As expected, the symbols on the maps corresponding to the MC-4 EeV model are mostly hexagons (Fe or sub-Fe nuclei), while many squares and pentagons (light and intermediate nuclei) contribute to the clustering on the maps corresponding to the MC-15 EeV model.

Going down to lower energies, the difference becomes obvious. In Fig. 21, we show typical sky maps for the same scenarios and the same realizations as above, with an energy threshold at 50 EeV. In the case of the MC-15 EeV model (right), many protons and light nuclei are observed near their sources (in red and in blue), while the Fe events are more loosely distributed in the case of the MC-4 EeV model (left). In passing, we note that the source whose events are colored in red in Fig. 19 appears in blue on the left sky map of Fig. 21, and vice-versa. Since the color code is associated with the rank of the source in terms of apparent luminosity, this means that the sources with the highest and the second highest luminosities have been interchanged between 50 EeV and 80 EeV. Although the sources are at different distances (35 Mpc and 44 Mpc), the difference between the effect of the GZK energy losses at these distances and energies is not significant. The main reason for the change in the relative luminosity is the magnification effect of the GMF, emphasized in Sect. 4.4.2. This effect is particularly strong for Fe nuclei in this energy range, owing to their intermediate rigidity. As a matter of fact, the source that appears in blue in Fig. 19 is located in a region where the magnification factor is greater than that of the red source, and even more so at rigidity R = 50 / 26 EV than at rigidity R = 80 / 26 EV. The source at 44 Mpc on the far left of the sky map is thus largely “boosted” at 50 EeV, compared to that at 35 Mpc in the middle right of the map. This is but one of the interesting features associated with the rigidity-dependent magnification effects (see also Sect. 4.4.2), whose consequences and astrophysical interest will be studied in a separate paper.

Coming back to the comparison between the MC-4 EeV and MC-15 EeV models, it is instructive to look at the fraction of the anisotropic realizations of these models as a function of the significance of their anisotropy (as measured by the two-point correlation function). This is shown in Fig. 22. The curves in the bottom panel correspond to the MC-15 EeV model and appear intermediate between those shown for the same model in Figs. 14 and 16 (right), as expected from the intermediate source density. The curves in the top panel correspond to the MC-4 EeV model and show that typically 5–10% fewer realizations display a given level of anisotropy with the JEM-EUSO statistics than in the previous case, and ~25% with the Pierre Auger Observatory statistics, owing to the absence of light nuclei. Overall, the same structure can nevertheless be observed and the same conclusions can be drawn.

However, in addition to the statistics and energy thresholds shown in this type of figure for the other models, we have shown in the two plots of Fig. 22 an additional curve corresponding to an extended Pierre Auger Observatory, which would collect 5000 events above 20 EeV (corresponding roughly to a total exposure of 100 000 km² sr yr). The difference between the MC-4 EeV and MC-15 EeV is striking at this energy, with these statistics. In the case of the MC-4 EeV, 5000 events above 20 EeV do not produce a more significant anisotropy than 84 events above 55 EeV or 67 EeV – or, more precisely, significant anisotropy is not obtained in a larger number of realizations of the MC-4 EeV model (with n_s = 10^-5 Mpc^-3) for 5000 events at 20 EeV than for 84 events at 55–67 EeV. In contrast, 5000 events above 20 EeV lead to an anisotropy as significant (or lead as often to a significant anisotropy) as the JEM-EUSO statistics above 50, 80, or 100 EeV, in the case of the MC-15 EeV model. This is due to the presence of protons in sufficient numbers at this low energy to produce small-scale anisotropies, despite the much more distant horizon and thus the larger number of sources contributing to the flux (see Lemoine & Waxman 2009).

This interpretation is confirmed by the curves plotted in Fig. 23, where we show the fraction of the realizations of the MC-4 EeV and MC-15 EeV models that display anisotropies with a significance of 2σ, 3σ, 4σ, with the extended Pierre Auger Observatory statistics at 20 EeV, as a function of angular scale. Clearly the angular scale where a large number of realizations of the MC-15 EeV model shows a large anisotropy corresponds to small deflections (thus associated with light nuclei), whereas this is not the case for the MC-4 EeV model, due to the dominance of heavier nuclei with larger deflections.

In conclusion, even though both scenarios are compatible with the current data, they could in principle be distinguished through their different anisotropy patterns as a function of energy with data sets collected by experiments having a larger exposure.

5.3. Proton-dominated models: MC-high and pure-p

In the previous section, we considered models in which the UHECR sources accelerate protons only up to an energy below the GZK energy scale. In these scenarios, the highest energy particles are heavy nuclei, experiencing large deflections as a result of their interactions with the magnetic field in our Galaxy. In contrast, if the sources are able to accelerate protons up to the highest energies observed, say up to 3 × 10²⁰ eV or above, then the UHECRs are likely to be dominated by protons in the GZK range. Indeed, the nuclei of hydrogen and helium are overwhelmingly dominant in the interstellar medium, so unless a strong discrimination mechanism favors heavier nuclei in the acceleration process, H and He nuclei should be dominantly injected in the intergalactic medium. Now, since the GZK horizon scale of the He nuclei (and of the light and intermediate nuclei) is much smaller than that of the protons in the energy range under consideration, the latter should dominate the UHECR composition observed on Earth.

The scenarios in which protons can be accelerated up to 300 EeV or above are thus particularly interesting in the context of the ongoing and very important quest for the UHECR sources. Figures 7, 8, and 10 clearly show that the deflections of protons in the GZK energy range are very small, except in a small part of the sky around the Galactic center. As a consequence, a proton-dominated scenario should easily lead to the detection of individual sources by direct pointing, once a few UHE protons are observed from the closest sources, at high enough energy for the roughly isotropic background from more distant sources to be negligible.

In Fig. 24, we show some typical sky maps (in the middle of the distribution, as defined above) obtained with the JEM-EUSO statistics at 50 EeV, 80 EeV, and 100 EeV, as well as with the current Pierre Auger Observatory statistics (top panel), for the proton-dominated model, MC-high, with a source density of 10^-4 Mpc^-3.

As expected, individual sources are easily detected in these scenarios. Unfortunately, although these sky maps appear to be a perfect situation for discovering and studying the UHECR sources, a simple look at the top panel of the figure shows that such a scenario is already excluded by the existing data. Indeed, the sky map built for the same realization with the Pierre Auger Observatory coverage and current statistics shows very tight multiplets, right in the direction of the sources. Such small-angular scale clusters of events have not been observed, whereas they are expected for 100% of the realizations of the MC-high scenario with n_s = 10^-4 Mpc^-3. The situation would, of course, be even worse if we drew the sky maps corresponding to an MC-high scenario with a lower source density, since the average multiplicities of the different sources would then be even higher and lead to even more obvious small-angular-scale multiplets that are impossible to miss in the Pierre Auger Observatory-like sky maps. As already indicated, the results obtained with a pure-proton model are exactly the same in all respects.

It is nevertheless interesting to see whether some specific scenarios in which protons dominate up to highest energies can be compatible with the Pierre Auger Observatory data. For this, we allowed the source density to be extremely high (see Pierre Auger Collaboration 2013) and simulated sky maps with a version of the MC-high model where all the galaxies in the catalog (see Sect. 3.3) were assumed to be sources of UHECRs. This corresponds to a density n_s = 1.6 × 10^-3 Mpc^-3.

A minority, but substantial fraction of the realizations of this scenario were indeed found to be compatible with the existing data, on the basis of the two-point correlation function. The corresponding probabilities are shown in the usual way in Fig. 25 (top). We also show examples of sky maps corresponding to this scenario in Fig. 26. The top panel shows the sky map of one of the realizations that gives one of the smallest anisotropy signals out of the 500 realizations simulated, with the reference Pierre Auger Observatory statistics and sky coverage. No obvious anisotropy can be seen on the map, which is indeed confirmed by the result of the statistical anisotropy study. From the point of view of the two-point correlation function, this sky map is similar to the actual Pierre Auger Observatory sky map at 55 EeV. The second panel of Fig. 26 shows the sky map obtained with the “median realization”, i.e., the realization sitting in the middle of the distribution of the 500 realizations ordered by significance of autocorrelation anisotropy. Finally, the third and fourth panels show the expected sky maps of this realization with the JEM-EUSO reference statistics with a threshold in the (reconstructed) energy of the UHECR events at 80 EeV and 100 EeV, respectively. These sky maps clearly show the interest of increasing the statistics for such scenarios. While the sources cannot be identified so far, they will certainly be so with a ten times greater exposure.

We performed simulations with an extended Pierre Auger Observatory statistics and found that significant multiplets should be detected by the Pierre Auger Observatory in the coming years as well. These results are summarized quantitatively in Fig. 25. As can be seen, 100% of the realizations give rise to extremely significant small-scale anisotropies at either the future Pierre Auger Observatory or JEM-EUSO statistics.

5.4. Pure-Fe models

We now turn to the study of an extreme case of heavy source composition, in which all the UHECRs injected in the intergalactic medium are Fe nuclei. As already noted, the GZK horizon structure for Fe nuclei is not very different from that of protons, so it is possible to obtain just as good a fit of the all-sky UHECR energy spectrum by using pure-Fe sources as by using a mixed-composition model, with either a high or a low value of the maximum proton energy (even though the required steepness of the source spectrum is different in each case, see above).

From what we know (or think we know) of the astrophysical environments where ultra-high-energy particle acceleration might occur in the universe, one must recognize, however, that a pure-Fe model is not very realistic by itself. Nevertheless, pure-Fe models, as a reference case of study, display interesting features that help in understanding some important aspects of the UHECR phenomenology, and are thus worth investigating.

This is mostly because UHE Fe nuclei, say around or above 10²⁰ eV, are photodissociated into lighter nuclei as they propagate through the intergalactic medium and interact with the background of microwave and infrared photons. These interactions typically eject one or a few nucleons out of the UHE nuclei, and since most of them occur near the threshold energy of the photodissociation (through giant dipolar resonance) in the rest frame of the nuclei, they roughly leave the remaining nuclei as well as the ejected nucleons with the same Lorentz factor, i.e., with the same energy per nucleon as the parent nuclei (see Allard 2012 for a more complete discussion of heavy nuclei propagation). As a consequence, an iron nucleus with initial energy E above the GZK cut-off will produce secondary particles and eventually protons with an energy E/A. By this process, the secondary protons coming from primary Fe nuclei can reach energies only a factor of 56 / 26 ≃ 2.15 lower than the energy whose primary protons could reach in the same sources if they were accelerated there, and if the maximum energy were proportional to the charge of the nuclei, Z, as assumed for the mixed-composition models.

As a consequence, pure-Fe models with a conceivably high maximum energy at the sources lead to UHECR populations that contain a non-negligible fraction of protons, when propagation effects are taken into account. For this reason, the pure-Fe models turn out to offer a simpler way to somehow explore complicated situations than the typical scenarios investigated here, where all the sources have the same spectrum, composition, and maximum energy.

In particular, we may use the pure-Fe models as a first hint of the phenomenology that could be relevant to scenarios in which, for example, the maximum energy of the protons is different in different sources, and can reach 10²⁰ eV (or more) in only a small number of sources (e.g., the most powerful ones). In such a case, most of the UHECRs in the GZK energy range would be heavy nuclei and thus experience relatively large deflections, as in the cases explored in Sect. 5.2, but protons and low-Z nuclei would also be present at some level, originating in a subset of the sources. Another possibility could be that the maximum energy reached in the UHECR sources is not as well-defined as usually assumed and does not give rise to an exponential cut-off. If the acceleration mechanism leads to a more gradual decrease in the UHECR flux injected in the extragalactic space as a function of energy, then a transition toward a heavy composition may occur, say, around 10¹⁹ eV (as possibly indicated by the Pierre Auger Observatory data) owing to the progressive extinction of the protons at the sources, while a subdominant fraction of protons survives up to the highest energies, with a decreasing fraction. Finally, a distribution of maximum energies has also been shown by Blaksley & Parizot (2012) to result in a modified apparent source composition, which may be important for understanding the evolution of the average UHECR mass as a function of energy. In such a situation, a transition toward heavier nuclei occurs at ultra-high energy, without implying a complete disappearance of protons or light/intermediate nuclei, as in the mixed-composition models with identical sources that are generally explored (e.g., MC-4 EeV and MC-15 EeV models).

For all these reasons, even though pure-Fe models are not astrophysically realistic by themselves, the results of this subsection should be considered as indicative of a number of interesting situations in which a small, but non-negligible fraction of protons are present among UHECRs up to the highest energies.

The distinctive feature of pure-Fe models (or more generally models with a subdominant fraction of protons, either primary or secondary) is evident on the sky maps shown in Fig. 27. While no significant clustering is observed in the data set corresponding to the current Pierre Auger Observatory statistics (84 events above 55 EeV, assuming the Pierre Auger Observatory energy scale), a number of UHECR sources appear very clearly with larger statistics, thanks to the secondary protons generated during propagation.

In the Pierre Auger Observatory reference sky map, two protons are seen in purple around their source (toward the top of the coverage map), but could be merely random coincidence. As for the most intense source in this sky map, located at 23 Mpc and marked by a red star at the bottom right of the map, it spreads many events across the sky, but only one of them is a proton. Although it is located very close to the source, it cannot be used to determine the source position, since it appears to be like any other events in the sky map.

It is worth noting, however, that the possibility of determining the composition of UHECRs event by event would be very important, especially for this kind of model, as a subsky map built by selecting only protons would give a view of the UHECR sky that is much easier to decipher. However, such a perspective appears quite distant at the moment, from the experimental point of view, be it only because of the so-called shower-to-shower fluctuations that can easily make a particular proton-induced atmospheric shower be very similar to an Fe-induced shower, and vice-versa. New detection techniques, e.g., based on the radio signal of the showers, may nevertheless change this situation in the future (Huege 2013; Revenu & Marin 2013).

The identification of proton events can, however, be done in an indirect, but rather trivial way without any detailed study of the showers, apart from the reconstruction of their arrival direction with a precision that is usually not experimentally challenging. Indeed, looking at the bottom panel of Fig. 27, one sees that the dominant source, which contributes 49 events above 100 EeV, spreads many UHE Fe nuclei over a large area of the map on the left hemisphere, as well as ten events in a relatively tight cluster on the right, but leads also to the superposition of nine events right in the direction of the source. It is statistically impossible for Fe nuclei to be observed so precisely in the same arrival direction, and since intermediate-mass nuclei (as can also be checked from the shape of the symbols on the map) are destroyed at lower energy by the GZK process, only protons can be responsible for such a tight cluster (barring a random coincidence, which may account for one or two events at most in a data set with such statistics).

The sky maps corresponding to an energy threshold of 50 EeV and 80 EeV (second and third panels of Fig. 27) also display some tight clusters, corresponding to the direction of individual sources. The evolution of the sky maps with energy can also give further information (and confidence) about the source location.

We have chosen here a realization that displays less significant anisotropies than the average, in order to remain compatible with the current constraints. We could also have presented a model with a higher source density, with essentially the same features. Indeed, many pure-Fe models already show significant anisotropies with the current Pierre Auger Observatory statistics, as is apparent in Fig. 28, where we show the usual statistical analysis of the anisotropy properties of the model. Respectively around 60% or 40% of the realizations appear to display small-scale anisotropies with a 2σ or 3σ significance, with the Pierre Auger Observatory statistics of reference.

6. Summary and perspective

In this paper, we have systematically studied some aspects of the phenomenology of a wide range of possible astrophysical scenarios for the origin of UHECRs, allowed by the current constraints and available data. We concentrated on the statistical analysis of the anisotropies that can be expected in the arrival directions of the particles at different energies, varying the source composition and spectrum, as well as the source density and the intergalactic magnetic field.

This study requires a proper treatment of the propagation of the UHECRs, whether protons or heavier nuclei, including the energy losses and photodissociation induced by the interaction with the various photon backgrounds, as well as the deflections caused by the interaction with the extragalactic and Galactic magnetic fields.

The main part of the angular distance between the source direction and the observed arrival direction of a given UHECR is caused by the interaction with the Galactic magnetic field, which is not fully characterized. We presented the model used for our computation and the resulting deflection properties as a function of rigidity in Sects. 4.3 and 4.4. Although the actual deflections are impossible to predict with confidence, we argued that the overall patterns relevant to our present study should be reliable on average, since they depend more on the typical structures present in the GMF than on the particular arrangement of the magnetic field maxima or reversals. This is also the reason we limited our anisotropy studies to a search for intrinsic anisotropies, which only depend on the global statistical properties of the GMF, rather than searching for features associated with particular directions in the sky. Our results are thus conservative in this respect, since other studies could in principle be carried out to derive additional constraints and information about the source models (e.g., by searching for multiplets with coherent patterns in the energy/deflection space or correlations with particular sources or structures).

Likewise, the exact distribution of the sources in the sky, as well as the UHECR deflections recorded as a consequence of the specific location of the sources with respect to the GMF structures cannot be considered as very important. By exploring many realizations of each scenario, we believe that the statistical properties that we obtained are robust to reasonable variations in the GMF model or other assumptions underlying our study. As discussed in Sect. 3.3, we used the 2MRS catalog as a reference catalog of possible sources, in order to respect the local distribution of matter as much as possible, in the absence of any information about the actual UHECR sources. But the main relevant parameter is the source density, and the exact location of the sources plays no role in the results presented here.

The astrophysical models we investigated were described and discussed in Sect. 3. They make use of earlier results, and explore situations in which the UHECRs are dominated by different types of particles (light, heavy, or with a transition from light to heavy), while ensuring compatibility with the existing constraints, notably as far as the UHECR energy spectrum and intrinsic anisotropies are concerned.

The general conclusion of our study is that the absence of a clear anisotropy signal in the current UHECR data does not imply that the UHECR sky will remain impossible to decipher and that UHECR sources cannot be found in a foreseeable future. On the contrary, we showed that a new generation of detectors, with a typical integrated exposure of 300 000 km² sr yr and full sky coverage, appears in a very good position to make a major breakthrough in this part of the field of astroparticle physics, by detecting significant anisotropies in essentially all the currently favored astrophysical scenarios.

The least favorable scenarios, in which protons are totally absent at the highest energies, have been shown to lead nevertheless to anisotropies with a significance of at least 4σ in the vast majority of the realizations. We would be able to derive important information about the UHECR sources with an instrument like JEM-EUSO or another large exposure detector except for very unlikely (and unlucky) distributions of the nearby sources. We made this statement quantitative by exploring the so-called cosmic variance for each astrophysical model through the simulation of 500 independent realizations of the different scenarios.

In more favorable scenarios, the presence of UHE protons can be extremely helpful for detecting the sources directly through identifying tight clusters of events (within a few degrees) in given directions of the sky. Although such scenarios have been challenged by the Pierre Auger Observatory results on UHECR composition, we also showed that proton-dominated scenarios are still viable from the point of view of anisotropies, provided that the source density is particularly large, on the order of 10^-3 Mpc^-3. For such scenarios, the value of a detector with larger statistics is particularly obvious, since the most nearby sources can be identified easily.

We pointed out that the generally studied scenarios, which lack more knowledge of the sources, intended to limit the number of free parameters, have a common drawback of assuming that all the sources have essentially the same parameters. As a matter of fact, we only studied standardized models, where all sources have the same composition and the same maximum energies. Although this allowed us to investigate the range of features that could be expected from a high-statistics, full-sky coverage UHECR detector in the future, the actual situation is most certainly more complicated, with some sources able to accelerate protons up to higher energies than others, so that a small, but crucial fraction of protons may still be present at high energy, even if the UHECRs turn out to be dominated by heavy nuclei, with smaller rigidities. A hint as to how this could modify our results has been given through study of pure-Fe models, which are characterized by a subdominant component of (secondary) protons. Qualitatively, it is clear that such more realistic models can only produce UHECR skies where more significant anisotropies can be detected and where individual sources can be identified with more precision, thanks to the possible existence of small angular scale multiplets, which are fully compatible with the current constraints derived from low-statistics data sets. Therefore, the situations considered here are in this respect conservative. A more exhaustive study of these more realistic astrophysical models, where different sources may have different values of E_max,p, is left for future work.

We also note that the tools presented in this paper can easily be applied to a wider variety of astrophysical scenarios, including those involving bursting sources, which have been left aside in the present study.

In this paper, we have used the JEM-EUSO performances as a reference for a future, large-exposure detector. This allowed us to be definite regarding the detection efficiency, which will most probably depend on energy for such large detectors. However, this choice is by no means limitative. Any other detector achieving the same exposure and performances (or better), would be equally (or more) suitable to write the next page of the long history of UHECR detection and cosmic ray studies. We made sure, however, that our results would properly apply to detectors with the somewhat lower performances that can be achieved in space (at least for a first generation of instruments). Indeed, a detector with a degraded precision on the reconstructed energy will attribute a given energy to a substantial number of events which are in reality less energetic. This so-called spill-over is asymmetric, owing to the steeply decreasing energy spectrum. Now, since lower energy events come from sources which are distributed over a larger volume than high-energy events, as a consequence of the GZK effect, a degraded energy resolution tends to degrade the anisotropy signal that can be detected at a given energy. To protect ourselves against a possible misinterpretation of our results as a consequence of this problem, we decided to apply to all the data sets generated in this study a (hopefully) conservative error on the energy reconstruction of 30%. In the case of a future detector with better energy resolution, the anisotropy signals to be expected can thus only be more significant than those presented here.

In conclusion, our results give strength to the general idea that a gain of about one order of magnitude in the total exposure of UHECR detectors would make a significant difference compared to what exists and allow considerable progress in the study of these mysterious particles, whose sources are still unknown and which challenge both the particle acceleration scenarios in the universe and the astrophysical modeling of the sources. The very likely detection of significant anisotropies (as made quantitative by our study) means that, for the first time, a situation could be reached in which it becomes useful (because meaningful) to study different regions of the sky separately, dominated by different UHECR sources. Further astrophysical information about the UHECR sources and acceleration processes could, for instance, be obtained by comparing the fluxes and spectra obtained in different regions of the sky and studying particular patterns in the UHECR energy/arrival direction space. These aspects of the problem have not been addressed here, since we focused on the first key question of anisotropies, but they should be investigated in future works.

Finally, while stressing the astrophysical value of pursuing UHECR studies with a new generation of detectors, we wish to call the importance of linking this study to the complementary studies of various classes of high-energy sources in the universe. UHECR studies need increased data sets, notably at the highest energies where the GZK effect ensures the dominance of only a handful of sources in the whole sky, but they will also greatly benefit from the development of high-energy astrophysics in general, with the central goal of applying multimessenger constraints to the study of what may still be regarded as the extreme universe.

Acknowledgments

The authors would like to thank Angela V. Olinto and Paolo Privitera for helpful discussions. We are grateful to Ronnie Jansson and Glennys R. Farrar for helping us implementing the Galactic magnetic field. We thank Nicolas G. Busca, Guillaume Decerprit, Noémie Globus and Kumiko Kotera for their involvement in previous works regarding extragalactic cosmic rays. Benjamin Rouillé d’Orfeuil was located at the Kavli Institute for Cosmological Physics at the University of Chicago when the project was initiated. His research was then supported by the National Science Foundation under grant No. NSF PHY-0551142 and a generous endowment from the Kavli Foundation.

References

All Tables

All Figures