format: text

text:
# Ultra-High Energy Cosmic Rays in a Structured and Magnetized Universe

## I Introduction

Over the last few years the detection of several giant air showers, either through ground based detectors haverah ; agasa  or fluorescence telescopes fe ; hires , has confirmed the arrival of ultra high energy cosmic-rays (UHECRs) with energies up to a few hundred EeV (1 EeV $`10^{18}`$eV). Their existence poses a serious challenge and is currently subject of much theoretical research as well as experimental efforts (for recent reviews see reviews ; bs-rev ; school ).

The problems encountered in trying to explain UHECRs in terms of “bottom-up” acceleration mechanisms have been well-documented in a number of studies (e.g., Refs. hillas-araa ; ssb ; norman ). In summary, apart from the specific energy draining interactions in the source the maximal UHECR energy is limited by the product of the accelerator size and the strength of the magnetic field. According to this criterion it turns out that it is very hard to accelerate protons and heavy nuclei up to the observed energies, even for the most powerful astrophysical objects such as radio galaxies and active galactic nuclei.

In addition, nucleons above $`70`$EeV suffer heavy energy losses due to photo-pion production on the cosmic microwave background (CMB) — the Greisen-Zatsepin-Kuzmin (GZK) effect gzk  — which limits the distance to possible sources to less than $`100`$Mpc stecker . Heavy nuclei at these energies are photo-disintegrated in the CMB within a few Mpc heavy . Unless the sources are strongly clustered in our local cosmic environment, a drop, often called the “GZK cut-off” in the spectrum above $`70`$EeV is therefore expected bbo , even if the injection spectra extend to much higher energies. However, the existence of the latter is not established yet from the observations mbo . In fact, whereas a cut-off seems consistent with the few events above $`10^{20}`$eV recorded by the fluorescence detector HiRes hires , it is not compatible with the 8 events (also above $`10^{20}`$eV) measured by the AGASA ground array agasa . The solution of this problem may have to await the completion of the Pierre Auger project auger  which will combine the two complementary detection techniques adopted by the aforementioned experiments.

Adding to the problem, there are no obvious astronomical counterparts to the detected UHECR events within $`100`$Mpc of the Earth elb-som ; ssb . At the same time, no significant large-scale anisotropy has been observed in UHECR arrival directions above $`10^{18}`$eV, whereas there are strong hints for small-scale clustering: The AGASA experiment has observed five doublets and one triplet within $`2.5^{}`$ out of a total of 57 events detected above 40 EeV agasa . When combined with three other ground array experiments, these numbers increase to at least eight doublets and two triplets within $`4^{}`$ uchihori . This clustering has a chance probability of less than $`1\%`$ in the case of an isotropic distribution.

Independent of the specific UHECR production mechanism, there are currently two possible explanations of the experimental findings described above: The first assumes very weak intergalactic magnetic fields capable of deflecting UHECRs only up to a few degrees, or neutral primaries. In this case the apparent isotropy would indicate that many sources contribute to the observed flux and most of these sources would be at cosmological distances because the local source distribution is in general too anisotropic to be consistent with the observed UHECR isotropy. This would also explain the absence of nearby counterparts and a subset of especially powerful sources would explain the small-scale clustering tt . Indeed, it has been argued that UHECR arrival directions correlate with the positions of BL Lacertae objects, suggesting these as sources accelerating protons bllac , although there seems to be disagreement about this in the literature bllac1 . Furthermore, some of these objects may be too far away to be consistent with the GZK effect, which would require new physics such as Lorentz symmetry violations jlm . In contrast, correlations with compact radio quasars have not been found c\_quasar . If correlations with astrophysical objects are confirmed, this would strongly suggest small deflection or neutral primary particles. Whatever the sources are in this scenario, for small deflection one can in principle constrain the characteristics of the magnetic fields along the line of sight and the source properties by analysing arrival times, directions, and energies of observed small-scale multi-plets sl . Also, in the small deflection scenario the experimental confirmation of a GZK cutoff is expected.

However, the assumption of weak intergalactic magnetic fields seems at odds with several observations kronberg . Most remarkable are the detections of Faraday rotation measures which seem to indicate field strengths at the $`\mu `$G level within the inner region ($``$ central Mpc) of galaxy clusters bo\_review . In addition, the recent mounting evidence for diffuse radio-synchrotron emission in numerous galaxy clusters gife00  and in a few cases of filaments kkgv89 ; bagchietal02 , seems to suggest the presence of magnetic fields as strong as 0.1-1.0$`\mu `$G at the relatively low density outskirts of collapsed cosmological structures. In fact, extragalactic magnetic fields (EGMF) as strong as $`1\mu G`$ in sheets and filaments of the large scale galaxy distribution, such as in our Local Supercluster, are compatible with existing upper limits on Faraday rotation bo\_review ; ryu ; blasi . It is also possible that fossil cocoons of former radio galaxies, so called radio ghosts, contribute to the isotropization of UHECR arrival directions mte . Thus, relatively strong magnetic fields seem to be ubiquitous in intergalactic space, although their theoretical understanding is still limited bt\_review .

Such observational evidence motivates a second, more realistic scenario, which takes into account the existence of strong ($`B0.11\mu `$G) intergalactic magnetic fields correlated with the large scale structure. In this case magnetic deflection of charged primaries would be considerable even at the highest energies and the observed UHECR flux could be dominated by relatively few sources within about 100 Mpc. Here, large scale isotropy could be explained by considerable angular deflection leading to diffusion up to almost the highest energies and the small scale clustering could be due to magnetic lensing hmrs . The locations of clusters of events of different energies would in this case coincide with the crossing points of the caustics for these energies where fluxes are enhanced.

In the present paper we take this second point of view and investigate in some detail the effects of propagation of UHECRs, assumed to be dominantly nucleons, in a magnetized large scale structure matter distribution computed according to a numerical cosmological simulation.

Early investigations of this scenario have been carried out in Refs. slb ; ils ; lsb ; sse ; is , assuming that sources and magnetic fields follow a pancake profile of scale height $`3`$Mpc and scale length $`20`$Mpc, the magnetic field having a power law spectrum at length scales below $`1`$Mpc. UHECR propagation was computed through a numerical code that accounts for magnetically induced deflections and all relevant energy losses slb ; lsb ; ils . The cases of a single source slb ; ils , as well as continuous lsb  and discrete source distributions is  have been investigated. The above studies led to the result that the multi-pole moments and autocorrelation functions of the arrival directions best fit the AGASA data for a number $``$10 sources in the Local Supercluster, assumed to emit continuously, and a maximal field strength of $`0.3\mu `$is .

Ideally, however, it would be desirable to study the propagation of UHECRs based on distributions of both potential sources and observed magnetic field properties. However, up to now, only catalogs of candidate sources have been available. Magnetic fields, on the other hand, have been approximated in a number of fashions: as negligible sommers , as uniform ynts , or as organized in spatial cells with a given coherence length and a strength depending as a power law on the local density tanco .

In the present paper we attempt to go beyond some of the above limitations by computing for the first time the propagation of the UHECRs in a magnetized cosmological environment computed through numerical simulations. We carry out a fully cosmological simulation of large scale structure formation which, in addition to dark matter and baryonic gas, follows the evolution of a passive magnetic field. This approach is motivated by the fact that $`\mu `$G magnetic fields are mostly negligible for the purpose of the dynamics of the large scale cosmic flows (hence their passive character). In addition, and basically for the same reason, the structure of magnetic fields on scales of interest for UHECR propagation ($``$ 100 kpc) is mostly determined by the hydrodynamic flow. This is confirmed by the fact that in these simulations, the magnetic field looses memory of its initial conditions, soon after the formation of structures begins. Finally, the statistical properties of cosmological structure in the universe are rather homogeneous. Therefore, the simulated matter structure and magnetic field distributions should provide a realistic scenario for studying the statistical properties of UHECR source distributions and propagation in a cosmic environment. In the present study we assume the sources to follow the baryon density. Furthermore, the observer is supposed to be in regions of the simulated matter distribution which contain structures of the same size and baryonic gas temperature as our local neighborhood. This should provide a suitable environment to simulate the arrival of UHECRs from extragalactic distances and the effects of local magnetic fields of various strengths.

In the future such studies can be further improved by computing constrained simulations that reproduce in detail the observed matter distribution of the local universe. Such a simulation has been used for the case of radio ghosts in Ref. mte  where, however, the magnetic fields were not followed but were rather assumed to scale with the gas density. Constrained simulations including magnetic fields are relevant for predicting quantitative features such as location of clustered events, phases of anisotropies etc. and will be used in a following study. We point out, however, that for the reasons given above, effects of the magnetic field and source distributions in the local universe should essentially be captured by the present approach at least up to “cosmic variance”. The latter represents variations due to different source and observer locations and will be estimated in our simulations.

We also restrict ourselves to UHECR nucleons, and we neglect the Galactic contribution to the deflection of UHECR nucleons since typical proton deflection angles in galactic magnetic fields of several $`\mu `$G are $`10^{}`$ above $`4\times 10^{19}`$eV medina , and thus in general are small compared to extra-galactic deflection in the scenarios studied in the present paper.

The simulation is described in more detail in the next section. There we also describe the general features of our method and define the statistical quantities used for comparison with the data. In Sect. 3 we present results and we conclude in Sect. 4.

## II Motivation and Outline of the Numerical Model

### II.1 Magnetic Deflection

Contrary to the case of electrons, for charged hadrons deflection is more important than synchrotron loss in the EGMF. To get an impression of typical deflection angles one can characterize the EGMF by its r.m.s. strength $`B`$ and a coherence length $`l_c`$. If we neglect energy loss processes for the moment, then the r.m.s. deflection angle over a distance $`rl_c`$ in such a field is $`\theta (E,r)(2rl_c/9)^{1/2}/r_L`$ wm , where the Larmor radius of a particle of charge $`Ze`$ and energy $`E`$ is $`r_LE/(ZeB)`$. In numbers this reads

$`\theta (E,r)`$ $``$ $`0.8^{}Z\left({\displaystyle \frac{E}{10^{20}\mathrm{eV}}}\right)^1\left({\displaystyle \frac{r}{10\mathrm{Mpc}}}\right)^{1/2}`$ (1)
$`\times \left({\displaystyle \frac{l_c}{1\mathrm{Mpc}}}\right)^{1/2}\left({\displaystyle \frac{B}{10^9\mathrm{G}}}\right),`$
for $`rl_c`$. This expression makes it immediately obvious why a magnetized Local Supercluster with fields of fractions of micro Gauss prevents a direct assignment of sources in the arrival directions of observed UHECRs; the deflection expected is many tens of degrees even at the highest energies. This goes along with a time delay

$`\tau (E,r)`$ $``$ $`r\theta (E,d)^2/4`$
$``$ $`1.5\times 10^3Z^2\left({\displaystyle \frac{E}{10^{20}\mathrm{eV}}}\right)^2\left({\displaystyle \frac{r}{10\mathrm{Mpc}}}\right)^2`$
$`\times \left({\displaystyle \frac{l_c}{\mathrm{Mpc}}}\right)\left({\displaystyle \frac{B}{10^9\mathrm{G}}}\right)^2\mathrm{yr},`$
which may be millions of years. A source visible in UHECRs today could therefore be optically invisible since many models involving, for example, active galaxies as UHECR accelerators, predict variability on much shorter time scales.

### II.2 Numerical Simulation of the Large Scale Structure

The formation and evolution of the large scale structure is computed by means of an Eulerian, grid based Total-Variation-Diminishing hydro+N-body code rokc93 . We adopt a canonical, flat $`\mathrm{\Lambda }`$CDM cosmological model with a total mass density $`\mathrm{\Omega }_m=0.3`$ and a vacuum energy density $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_m=0.7`$. We assume a normalized Hubble constant $`h_{67}H_0/67`$ km s-1 Mpc-1 = 1 and a baryonic mass density, $`\mathrm{\Omega }_b=0.04`$. The simulation is started at redshift $`z60`$ with initial density perturbations generated as a Gaussian random field and characterized by a power spectrum with a spectral index $`n_s=1`$ and “cluster-normalization” $`\sigma _8=0.9`$.

We adopt a computational box size of $`50h_{67}^1`$Mpc. In this box the dark matter component is described by 2563 particles whereas the gas component is evolved on a comoving grid of 5123 zones. Thus each numerical cell measures about $`100h_{67}^1`$kpc (comoving) and each dark matter particle corresponds to $`2\times 10^9h_{67}^1M_{}`$. Besides the box and dark matter particle sizes the cosmological simulation is the same as that presented in Ref. miniati .

The magnetic field is followed as a passive quantity, that is magnetic forces are neglected. This is consistent with the strength of observed magnetic fields in most diffuse extragalactic environments. Basically we solve the induction equation with the velocity field provided by the simulated flow kcor97  and the initial magnetic field seeds generated by the Biermann battery mechanism. However, as already pointed out, the initial conditions are not important as the topological properties of the magnetic field are determined by the subsequent evolution of the large scale flow. This is responsible for its amplification through gas compression and shear flows. Thus, at the end of the simulation, the relative strength of the magnetic field in different regions is determined by the hydrodynamic properties of the flow. While the simulation outcome regarding the relative magnetic field strength and topology distribution are obviously retained, the overall normalization is chosen in order to reproduce the fields of several micro Gauss observed in the regions of largest density, namely galaxy cluster cores. Fig. 1 illustrates an example of the simulated magnetic pressure (top) and baryonic density (bottom) distributions. The figure shows two-dimensional cuts corresponding to a depth of $`100h_{67}^1`$kpc. The color images are in log scale and, for visualization purposes, span a dynamic range of 3 and 6 orders of magnitude for magnetic pressure and baryonic density respectively. The magnetic field is particularly strong in both postshock regions and inside relatively large structures where it has been compressed and stretched. Apparently, its distribution is less concentrated than the baryonic density, resembling in this respect that of the thermal pressure (not shown).

### II.3 Simulated UHECR Experiments

To simulate the propagation and arrival of UHECRs in the computational box we need to choose: (a) the location of the observer and (b) the source distribution. As anticipated in the introduction, the location of the observer is identified as a region whose general features in terms of scale, mass and temperature, resemble those of the local universe. That means a small group of galaxies characterized by a gas temperature of order of a fraction of a keV. There are several such structures in a $`50h_{67}^1`$Mpc box such as the one employed here. In the neighborhood of the one we selected as the observer location, we also find a larger group of galaxies with temperature of a few keV. In order to orient the simulation box with respect to the observed sky, the latter object, located at a distance of $``$ 34 Mpc, is arbitrarily associated with the Virgo cluster. This reference frame allows us to define a celestial system of coordinates $`(\alpha ,\delta )`$ which describes the arrival direction of events recorded by our virtual observer. It will be useful in the next section where the arrival direction probability distribution is constructed. The above setting is sufficient for the current purpose of investigating the effects on the propagation of UHECRs of realistic, topologically structured magnetic fields of various strengths.

We then chose randomly a certain total number $`N_s`$ of sources in the box, corresponding to an average source density $`8\times 10^6h_{67}^3N_s\mathrm{Mpc}^3`$, with probability proportional to the local baryon density. In order to avoid introducing too many free parameters, we further assume that all sources roughly emit the same power law spectrum of CRs extending up to $`10^{21}`$eV, with roughly equal total power. We also assume that neither total power nor the power law spectral index change significantly on the time scale of UHECR propagation. This can be up to a few Giga years for the magnetic fields considered here. Injected power and spectral index are then treated as parameters which can be fit to reproduce the observed spectrum, as will be seen below.

For each such configuration many nucleon trajectories originating from the sources were computed numerically by solving the equation of motion for the Lorentz force and checking for pion production every fraction of a Mpc according to the total interaction rate with the CMB and, in case of an interaction, by randomly selecting the secondary energies according to the differential cross section. Pair production by protons is treated as a continuous energy loss process.

A detection event was registered and its arrival directions and energies recorded each time the trajectory of the propagating particle crossed a sphere of radius 1 Mpc around the observer. For each configuration this was done until 5000 events where registered. For more details on this method see Refs. slb ; lsb ; ils .

### II.4 Data Processing

For each realization of sources and observer, these events were used to construct arrival direction probability distributions, taking into account the solid-angle dependent exposure function for the respective experiment and folding over the angular resolution.

For the exposure function $`\omega (\delta )`$ we use the parameterization of Ref. sommers  which depends only on declination $`\delta `$,

$`\omega (\delta )`$ $``$ $`\mathrm{cos}(a_0)\mathrm{cos}(\delta )\mathrm{sin}(\alpha _m)+\alpha _m\mathrm{sin}(a_0)\mathrm{sin}(\delta ),`$
$`\text{where}\text{ }\alpha _m=\{\begin{array}{cc}0\hfill & \text{if }\xi >1\hfill \\ \pi \hfill & \text{if }\xi <1\hfill \\ \mathrm{cos}^1(\xi )\hfill & \text{otherwise}\hfill \end{array},`$
$`\text{with}\text{ }\xi {\displaystyle \frac{\mathrm{cos}(\theta _m)\mathrm{sin}(a_0)\mathrm{sin}(\delta )}{\mathrm{cos}(a_0)\mathrm{cos}(\delta )}}.`$
For the AGASA experiment $`a_0=35^{}`$, $`\theta _m=60^{}`$, and the angular resolution $`2.4^{}`$ are used. For a full-sky Pierre Auger type experiment we add the exposures for the Southern Auger site with $`a_0=35^{}`$ and a putative similar Northern site with $`a_0=39^{}`$, and $`\theta _m=60^{}`$ in both cases, with an assumed angular resolution of $`1^{}`$.

From the distributions obtained in this way typically 1000 mock data sets consisting of $`N`$ observed events were selected randomly. For each such mock data set or for the real data set we then obtained estimators for the spherical harmonic coefficients $`C(l)`$ and the autocorrelation function $`N(\theta )`$. The estimator for $`C(l)`$ is defined as

$$C(l)=\frac{1}{2l+1}\frac{1}{𝒩^2}\underset{m=l}{\overset{l}{}}\left(\underset{i=1}{\overset{N}{}}\frac{1}{\omega _i}Y_{lm}(u^i)\right)^2,$$
(7)
where $`\omega _i`$ is the total experimental exposure at arrival direction $`u^i`$, $`𝒩=_{i=1}^N1/\omega _i`$ is the sum of the weights $`1/\omega _i`$, and $`Y_{lm}(u^i)`$ is the real-valued spherical harmonics function taken at direction $`u^i`$. The estimator for $`N(\theta )`$ is defined as

$$N(\theta )=\frac{C}{S(\theta )}\underset{ji}{}\left\{\begin{array}{cc}1\hfill & \text{if }\theta _{ij}\text{ is in same bin as }\theta \hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right\},$$
(8)
and $`S(\theta )`$ is the solid angle size of the corresponding bin. In Eq. (8) the normalization factor $`C=\mathrm{\Omega }_e/(N(N1))`$, with $`\mathrm{\Omega }_e`$ denoting the solid angle of the sky region where the experiment has non-vanishing exposure, is chosen such that an isotropic distribution corresponds to $`N(\theta )=1`$.

The different mock data sets in the various realizations yield the statistical distributions of $`C(l)`$ and $`N(\theta )`$. One defines the average over all mock data sets and realizations as well as two errors. The smaller error (shown to the left of the average in the figures below) is the statistical error, i.e. the fluctuations due to the finite number $`N`$ of observed events, averaged over all realizations. The larger error (shown to the right of the average in the figures below) is the “total error”, i.e. the statistical error plus the cosmic variance. Thus, the latter includes the fluctuations due to finite number of events and the variation between different realizations of observer and source positions.

Given a set of observed and simulated events, after extracting some useful statistical quantities $`S_i`$, namely $`C_l`$ and $`N(\theta )`$ defined above, we define

$$\chi _n\underset{i}{}\left(\frac{S_{i,\mathrm{data}}\overline{S}_{i,\mathrm{simu}}}{\mathrm{\Delta }S_{i,\mathrm{simu}}}\right)^n,$$
(9)
where $`S_{i,\mathrm{data}}`$ refers to $`S_i`$ obtained from the real data, and $`\overline{S}_{i,\mathrm{simu}}`$ and $`\mathrm{\Delta }S_{i,\mathrm{simu}}`$ are the average and standard deviations of the simulated data sets. This measure of deviation from the average prediction can be used to obtain an overall likelihood for the consistency of a given theoretical model with an observed data set by counting the fraction of simulated data sets with $`\chi _n`$ larger than the one for the real data.

## III Results

In the following we compare the results obtained for the simulated UHECR propagation experiments described above with the observational results. In accord with what was outlined in the previous section, the comparison is based on the statistical properties of the simulated and observed events, expressed in terms of the angular power spectrum and the autocorrelation function of the UHECR arrival distributions. A summary of the simulations run is contained in Tab. 1. There, for comparison, simulations 2 and 5 were performed for an observer situated in a small void with weak ambient magnetic fields.

We find that as long as the observer is surrounded by magnetic fields of about $`0.1\mu `$G, $`N_s10`$ nearby sources, i.e. sources within the simulation box, are necessary to reproduce multi-poles and autocorrelations marginally consistent with present data, limited, we emphasize, to the Northern hemisphere only. However, consistency of large scale multi-poles is somewhat worse than for the spatially more extended EGMF assumed in previous work is . In Figs. 2 and 3 we show as an example the results for the case of $`N_s=100`$ nearby sources, scenario 1 in Tab. 1, corresponding to a source density of $`8\times 10^4h_{67}^3\mathrm{Mpc}^3`$. The overall likelihood for $`n=4`$ in Eq. (9) is $`0.13`$ and $`0.63`$ for the multi-poles and autocorrelations shown, respectively. Also Fig. 4 shows that, for UHECR sources characterized by a proton injection spectrum roughly as $`E^{2.4}`$ and extending up to $`10^{21}`$eV, the observed spectrum at sub-GZK energies is well reproduced. In addition, above GZK energies the spectral slope is predicted to be somewhere between the AGASA and HiRes observations, see Fig. 4. Normalizing to the observed flux results in a UHECR power of $`5\times 10^{41}\mathrm{erg}\mathrm{s}^1`$ per source to be continuously emitted above $`10^{19}`$eV.

The situation for $`N_s=10`$ nearby sources does not lead to significantly different likelihoods, see scenarios 3 and 4 in Tab. 1. However, the case of just one source is clearly disfavored in terms of the multi-poles, see scenario 6 in Tab. 1. This confirms similar findings in earlier work ils .

If the observer is in a region of EGMF strength much smaller than $`0.1\mu `$G, as in scenario 2 of Tab. 1, for $`N_s100`$ nearby sources the predicted UHECR sky distribution reflects the highly structured large scale galaxy distribution, smeared out only by the fields surrounding the sources. This becomes obvious from Fig. 5 which shows that UHECR arrival directions are much less isotropic in this case than if the observer is immersed in fields $`B0.1\mu `$G.

Nevertheless, the overall likelihood significance for multi-poles up to $`l=10`$ is $`0.1`$, and thus not significantly worse than for the strong observer field case of Fig. 2. Therefore, the number of events observed by AGASA above $`40`$EeV is insufficient to distinguish this low observer field case from the strong observer field case based on anisotropy alone. However, as can be seen from Fig. 6, the low observer field case results in auto-correlations at angles $`\theta 3^{}`$ much larger than observed by AGASA. This is because strong magnetic fields at the observer position cause enough UHECR diffusion that their large-scale auto-correlations are significantly suppressed, as in Fig. 3. However, for fields considerably larger than $`0.1\mu `$G the auto-correlations tend to become too strong again, see scenario 4 in Tab. 1, probably due to increased magnetic lensing.

We also find that sources outside our Local Supercluster do not contribute significantly to the observable flux if the observer is immersed in magnetic fields above about $`0.1\mu `$G and if the sources reside in magnetized clusters and super-clusters: For particles above the GZK cutoff this is because sources outside the Local Supercluster are beyond the GZK distance. On the other hand, sub-GZK particles are mainly confined in their local magnetized environment and thus exhibit a much higher local over-density than their sources. Further, the suppressed flux of low energy particles leaving their environment is largely kept away from the observer if he is surrounded by significant magnetic fields is . Both effects can be understood qualitatively by matching the flux $`j(E)`$ in the unmagnetized region with the diffusive flux $`D(E)n(E,𝐫)`$ in terms of the diffusion coefficient $`D(E)`$ and the density $`n(E,𝐫)`$ of particles of energy $`E`$ which shows that the density gradient always points to the source. More quantitatively, the shape of the large-distance component is demonstrated in Fig. 7 which shows the observable flux resulting from an $`E^{2.4}`$ spectrum injected isotropically at a sphere with a radius of 40 Mpc around the observer. Note that despite the smaller energy losses the sub-GZK particles arriving from outside the Local Supercluster are likely to have a spectrum even more strongly suppressed than in Fig. 7 at low energies due to their containment in the source region. A significant contribution from sources at cosmological distances to sub-GZK energies thus requires that neither these sources nor the observer are immersed in too strong magnetic fields and/or an injection spectrum considerably steeper than $`E^{2.4}`$ to compensate for the systematic suppression of flux of lower energy particles.

The confidence levels that can be obtained with this method for specific models of our local magnetic and UHECR source neighborhood will greatly increase with the increase of data from future experiments. Full sky coverage alone will play an important role in this context as many scenarios predict large dipoles for the UHECR distribution. This is the case for basically all scenarios considered here, as demonstrated in Fig. 8. Whereas current northern hemisphere data are consistent with scenarios with $`N_s10`$ nearby sources at the $`1.5`$sigma level if the observer is surrounded by relatively strong fields $`B0.1\mu `$G, a comparable or larger exposure in the southern hemisphere would be sufficient in these cases to find a dipole at several sigma confidence level, as demonstrated in Fig. 8.

Finally, the distributions of events down to $`10^{19}`$eV also contain important information. Fig. 9 shows the multi-poles predicted by our standard scenario 1 in Tab. 1 that full-sky experiment would observe for 1500 events detected above $`10^{19}`$eV. This corresponds to twice the number of currently observed AGASA events and thus approximately reflects the current exposure. A corresponding figure for the AGASA detector alone would look similar. It is obvious that there is significant anisotropy even at $`l10`$, inconsistent with current AGASA observations. On the other hand, cosmic variance becomes more important at these lower energies, and a possible significant contribution from large-distance sources cannot be excluded if their magnetization is not too high, as discussed above. It is easy to see from Eq. (7) that if a fraction $`f_a`$ of $`N`$ events observed stems from an anisotropic, local contribution, whereas the fraction $`1f_a`$ is cosmological and completely isotropic, then

$$C_lC_{l,i}\left((1f_a)^2+\frac{C_{l,a}}{C_{l,i}}f_a^2\right),$$
(10)
where $`C_{l,i}=(4\pi N)^1`$ and $`C_{l,a}`$ are the expectation values of $`C_l`$ for the isotropic and the anisotropic distribution, respectively. Therefore, at $`10^{19}`$eV an isotropic cosmological flux about a factor 3 higher than the anisotropic flux originating within $`50`$Mpc would be needed to explain the isotropy observed by AGASA. For charged primaries this implies steep injection spectra and/or weak magnetic fields around observer and sources, as explained above. Without going into a more detailed analysis we remark that this will also require to decrease the flux contribution from nearby sources shown in Fig. 4 at the low energy end. As a consequence, the best fit injection spectrum for the local component will be slightly harder than the power law indices shown in Tab. 1. This is consistent with what is expected from shock acceleration theory rel\_shock .

## IV Conclusions

In the present work we performed UHECR propagation simulations based on the distributions of magnetic field and baryon density obtained from a simulation of large scale structure formation. The magnetic field was simulated as a passive quantity and normalized at simulation end in agreement with published measurements of Faraday rotation measures for groups and clusters of galaxies bo\_review . We considered finite numbers of discrete UHECR sources with equal total power and injection spectrum. Their positions were randomly selected with probability proportional to the baryon density. The observer was chosen within small groups of galaxies characterized by gas temperatures around a fraction of a keV, typical for our local environment. One chosen observer was found in a relatively high field region with $`B0.1\mu `$G. For comparison, we also chose an observer situated in a small void, where the surrounding field is $`B10^{11}`$G. We found that good fits to the AGASA data above $`4\times 10^{19}`$eV in the Northern hemisphere are only obtained for $`N_s10`$ sources and for observers surrounded by $`0.1\mu `$G fields. Otherwise the predicted arrival direction distribution is either too anisotropic or produces too large auto-correlations at angles larger than a few degrees. The best fit case occurs for $`N_s100`$, significantly higher than in previous work is  due to the more localized and more strongly structured magnetic fields considered here.

For the required local source number density and continuous power per source we find $`n_{\mathrm{source}}10^410^3h_{67}^3\mathrm{Mpc}^3`$, and $`Q_{\mathrm{source}}5\times 10^{41}\mathrm{erg}\mathrm{s}^1`$ respectively, the latter within about one order of magnitude uncertainty to both sides. This corresponds to an average UHECR emissivity of $`q_{\mathrm{UHECR}}=n_{\mathrm{source}}Q_{\mathrm{source}}10^{38}\mathrm{erg}\mathrm{Mpc}^3\mathrm{s}^1`$ also with an uncertainty of roughly one order of magnitude, not larger, since it is fixed by the observed UHECR flux.

Possible sources marginally consistent with these energy requirements are radio galaxies. Their present energy release of $`q_{\mathrm{rg}}5\times 10^{39}\mathrm{erg}\mathrm{s}^1\mathrm{Mpc}^3(f_{\mathrm{power}}/10)`$ ebkw97  is roughly what is required in order to produce a sufficient flux of UHECR, assuming that the injection power law is flat ($`E^2`$) rb93 ; rsb93 . The parameter $`f_{\mathrm{power}}`$ describes the ratio of the total power of the radio galaxy to the equipartition estimate based on its radio luminosity, and it enters the used radio luminosity-jet-power relation of ebkw97 . We expect $`f_{\mathrm{power}}10`$ within an order of magnitude. In the estimate of the radio galaxy power the observed radio luminosity function dp90  was integrated only for sources with a 2.7 GHz luminosity of more than $`L_{\mathrm{min}}=2\times 10^{22}\mathrm{Watt}\mathrm{Hz}^1`$, since they correspond to a luminosity of $`10^{43}\mathrm{erg}\mathrm{s}^1(f_{\mathrm{power}}/10)`$. This would provide the required UHECR luminosity per source of $`Q_{\mathrm{source}}5\times 10^{41}\mathrm{erg}\mathrm{s}^1(f_{\mathrm{power}}f_{\mathrm{UHECR}}/0.5)`$, using the optimistic assumption of rb93 ; rsb93  that $`f_{\mathrm{UHECR}}=310\%`$ of the radio galaxy power is converted into UHECRs. The implied number density of these radio galaxies is $`n_{\mathrm{source}}10^4\mathrm{Mpc}^3`$ and, therefore, is only marginally consistent with the required $`n_{\mathrm{source}}10^410^3h_{67}^3\mathrm{Mpc}^3`$. Since the number density increases strongly with decreasing $`L_{\mathrm{min}}`$ this requirement can possibly be fulfilled by allowing for a larger number of less powerful UHECR sources. This implies that basically every radio galaxy has to be an efficient UHECR source, not only the most powerful ones. Since many of the weaker radio galaxies do not exhibit a hot spot, which is assumed to be the UHECR acceleration site in the scenario of rb93 ; rsb93 , their efficiency in producing UHECR might be largely reduced. This is a potential serious problem for this scenario, since lowering $`f_{\mathrm{UHECR}}`$ by several orders of magnitude can not be fully compensated by assuming a higher radio galaxy jet-power, because $`f_{\mathrm{power}}>100`$ does not seem to be consistent with observations of radio galaxies 1992A&A…265….9F .

To conclude, radio galaxies can be the sources of UHECRs if even weak radio galaxies are efficient particle accelerators to the highest energies, otherwise they have serious problems to reproduce the smooth UHECR arrival direction distribution.

We also found that consistency with the isotropy observed by AGASA down to $`10^{19}`$eV requires the existence of an isotropic component with a flux about a factor 3 larger than the local component. This isotropic component would presumably be of cosmological origin and thus would not contribute significantly above $`4\times 10^{19}`$eV due to the GZK effect, consistent with the fact that at these energies we find local scenarios consistent with all data. The resulting best fit injection spectrum for the local component is $`E^{(2.22.4)}`$. In contrast, for the charged primaries of the cosmological component to dominate around $`10^{19}`$eV steep injection spectra and/or weak magnetic fields around observer and sources would be required. These two conflicting requirements provide a strong argument against the hypothesis of local astrophysical sources of UHECRs above $`10^{19}`$eV in a strongly magnetized and structured intergalactic medium.

Finally, we have also demonstrated that already a modest increase in data together with full-sky coverage will allow to put considerably stronger constraints on UHECR source and magnetic field scenarios than presently possible. In particular, our local scenarios predict the emergence of significant dipoles and quadrupoles above $`4\times 10^{19}`$eV.

Modeling our cosmic neighborhood and simulating UHECR propagation in this environment will therefore become more and more important in the coming years. This will also have to include the effects of the Galactic magnetic field and an extension to a possible heavy component of nuclei. For first steps in this direction see, e.g. Refs. mte ; ames , and Ref. bils , respectively.

## Acknowledgments

We would like to thank Martin Lemoine and Claudia Isola for earlier collaborations on the codes partly used in this work and to F. W. Stecker for useful comments to the manuscript. The work by FM was partially supported by the Research and Training Network “The Physics of the Intergalactic Medium” set up by the European Community under the contract HPRN-CT2000-00126 RG29185. The computational work was carried out at the Rechenzentrum in Garching operated by the Institut für Plasma Physics and the Max-Planck Gesellschaft.