text:
# Small Scale Clustering in the Isotropic Arrival Distribution of Ultra-High Energy Cosmic Rays and Implications for Their Source Candidates

## 1 INTRODUCTION

In spite of a increasing amount of data, the origin of cosmic rays in particular at the highest energy ($`10^{20}`$eV) is still unknown. It is very difficult to produce such ultra-high energy cosmic rays (hereafter UHECRs) in astrophysical shocks, which are thought to be responsible for the galactic cosmic rays. The gyroradii of UHECRs with energies of 1020eV are larger than our Galaxy, and this suggests an extragalactic origin for any astrophysical scenario. Using a Hillas-plot (Hillas, 1984; Selvon, 2000), active galactic nuclei (AGN), gamma-ray bursts (GRBs), dead quasars and colliding galaxies are considered as probable candidates. These conventional acceleration models are called bottom-up scenarios. On the other hand, this situation has in recent years triggered a number of exotic production models based on new physics beyond the Standard model of the particle physics (see Bhattacharjee and Sigl, 2000, and references therein), which are called top-down scenarios. In such scenarios, UHECRs owe their origin to decay of some supermassive particles, which could be produced from Topological Defects (TDs), or be certain metastable supermassive relic particles (MSRPs) produced in the early universe.

In bottom-up scenarios in which protons are accelerated in the potential UHECRs sources, like GRBs and AGNs, the large distances to the earth lead to another problem, called the GZK effect (Greisen, 1966; Zatsepin and Kuz’min, 1966). UHE protons with energy above 8$`\times 10^{19}`$ eV interact with photons of Cosmic Microwave Background (CMB) and lose their energy rapidly due to the production of pions over distances of tens of Mpc (Stanev et al., 2000). This should results in a cut-off of the cosmic ray spectrum at energy around 8$`\times 10^{19}`$ eV. The situation for nuclei is considered to be worse due to the photo-disintegration mechanism (Stanev et al., 2000). Observationally, there seems to be a disagreement between the AGASA which revealed the extension of the cosmic ray energy spectrum above $`10^{20}`$ eV (Takeda et al., 1998) and the High Resolution Fly’s Eye (HiRes; Wilkinson et al., 1999) which recently reported the cosmic ray flux with the GZK cut-off (Abu-Zayyad et al., 2002).

Since UHECRs above $`10^{20}`$ eV observed by the AGASA must originate within the GZK sphere, whose radius is typically about 50 Mpc, the arrival directions are expected to point toward their sources within a few degrees, which is typical deflection angle of UHECRs (e.g., Blasi and Olinto, 1998). However, no plausible source counterparts have been found within the GZK sphere. There is also a problem concerning the arrival directions of UHECRs (e.g., Blasi and Olinto, 1998). Their distribution seems to be isotropic on a large scale with a statistically significant small scale clustering (Takeda et al., 1999). The current AGASA data set of 57 events above 4 $`\times 10^{19}`$ eV contains four doublets and one triplet within a separation angle of 2.5. Chance probability to observe such clusters under an isotropic distribution is about 1 $`\%`$ (Hayashida et al., 2000). Potential models of UHECR origin are constrained by their ability to reproduce the observed energy spectrum and the isotropic arrival distribution of UHECRs with the small scale clustering.

A number of studies has been devoted to the problem of the energy spectrum. Yoshida and Teshima (1993) show the GZK cutoff in consideration of evolution of the universe. Blanton et al.(2001) show that the observed number of events above $`10^{20}`$ eV is 2 $`\sigma `$ higher than the expected one, assuming a hard $`(E^2)`$ injection spectrum and small local overdensity (a factor of two). Berezinsky et al.(2002) analyze three source models, GRB, AGN and local enhancement of UHECR sources. They conclude the AGASA events above $`10^{20}`$ eV, if confirmed, must be of different origin. Although these authors show that it is difficult to reproduce the observed extension of the energy spectrum by the bottom-up scenarios, they do not take into account of the effects of the extragalactic magnetic field (EGMF) on the energy spectrum. Furthermore, since they do not consider angular source distribution as well as the EGMF, the problems concerning the arrival distribution of UHECRs can not be discussed.

Although arrival distributions from a single source are calculated in the presence of magnetic field in some papers (Lemoine et al., 1997; Stanev et al., 2000), it appears to be difficult to explain the observational isotropy by a single source. With consideration of both the magnetic field and the source distribution, the energy spectrum and the distribution of the arrival directions are calculated in Sigl et al.(1999) and Lemoine et al.(1999). However the source distribution seems to be too simple, a two dimensional sheet (Sigl et al., 1999) or a pancake-like profile (Lemoine et al., 1999), which represent the Local Super Cluster (LSC).

Smialkowski et al.(2002) calculate arrival distribution of UHECRs using the IRAS PSCz catalog of IR galaxies as a realistic source model. However, they restrict their attention to galaxies with far infrared luminosity $`L_{\mathrm{FIR}}>10^{11}L_{}`$ triggered by collision and merging processes, which are possible sites of UHECR acceleration. It is also suspected that this galaxy sample do not contain some nearby galaxies in the LSC. Since nearby sources within the GZK sphere are responsible for UHECRs above $`10^{20}`$ eV, galaxy sample with better completeness on nearby galaxies should be used.

Given these situations, we examined in our previous paper (Ide et al., 2001) whether the current observation can be explained by a bottom-up scenario in which source distribution of UHECRs is proportional to that of galaxies. We used the Optical Redshift Survey (ORS; Santiago et al., 1995) as a realistic source model. Completeness on nearby galaxies of this sample would be better than that of the IRAS PSCz. By using the ORS sample, we can accurately quantify the contribution of nearby sources to the energy spectrum above $`10^{20}`$ eV. We also took effects of the EGMF on the energy spectrum and the arrival distribution of UHECRs into account.

However, there are problems in our previous study. First, we did not correct absence of galaxies in the zone of avoidance ($`|b|<20^{}`$) when using the ORS sample, so that we could not present accurate discussion of statistics on the arrival directions of UHECRs. Second, we considered propagation of UHECRs only over 40 Mpc, focusing on UHECRs with energies above $`8\times 10^{19}`$ eV. In this case, contribution of distant sources (outside 40 Mpc) to the cosmic ray flux below $`8\times 10^{19}`$ eV can not be taken into account. Finally, we injected UHECRs from all the galaxies including dwarf ones. Luminous galaxies are generally larger than dwarf ones, and thus they are expected to be able to accelerate particles to higher energy. Thus, we should investigate dependence of the results on the limiting magnitude of galaxies.

This paper is an extended study of our earlier work. We calculate the energy spectrum and the arrival distribution of UHECRs using the ORS galaxy sample. We correct the sample taking the selection effect and absence of galaxies in the zone of avoidance ($`|b|<20^{}`$) into account. This allow us to accurately discuss statistics on the arrival directions of UHECRs. We also extend the propagation distance of UHECRs, from 40 Mpc to 1 Gpc. As s result, we can calculate the cosmic ray flux below $`8\times 10^{19}`$ eV, properly including contribution from distant sources. Furthermore, we investigate dependence of the results on the limiting magnitude of galaxies.

We compare the results of the energy spectra with the current observations, AGASA and HiRes. In order to explore the possibility of reproducing the extension of energy spectrum beyond $`10^{20}`$ eV, we adopt several strength of the EGMF (1,10 and 100 nG). We also investigate dependence of the arrival distribution of UHECRs on the limiting magnitude of galaxies and the source number, and estimate the most favorable ones from comparison with the AGASA data using the harmonic analysis and the two point correlation function as statistical methods. We show that the three observable quantities including the GZK cutoff of the energy spectrum can be reproduced in the case that the number fraction $`10^{1.7}`$ of the ORS galaxies more luminous than $`20.5`$ mag is selected as UHECR sources. Implications of the results for the UHECR source candidate are discussed in detail.

In section 2, we show our method of calculation. Results are shown in section 3 and their implications are discussed in section 4. We conclude in section 5.

## 2 METHOD OF CALCULATION

### 2.1 Method of Calculation for Propagation of UHECRs

This subsection provides the method of Monte Carlo simulations for propagating protons in intergalactic space. At first, we assume that the composition of UHECRs is proton, and inject an $`E^2`$ spectrum within the range of ($`10^{19.5}`$ \- $`10^{22}`$)eV. 10000 protons are injected in each of 26 energy bins, that is, 10 bins per decade of energy for the cases of $`l_\mathrm{c}`$ = 40Mpc and ($`B`$, $`l_\mathrm{c}`$) = (1,10nG, 10Mpc), where $`B`$ and $`l_\mathrm{c}`$ are strength and correlation length of the extra galactic magnetic field (explained below). For another cases, we inject 5000 protons in each of the bins.

Particles below $`8\times 10^{19}`$ eV lose their energies mainly by pair creations and above it by photopion production (Yoshida and Teshima, 1993) in collisions with photons of the CMB. The pair production can be treated as a continuous loss process considering its small inelasticity ($`10^3`$). We adopt the analytical fit functions given by Chododowski et al. (1992) to calculate the energy loss rate for the pair production on isotropic photons. According to them, the energy loss rate of a relativistic nucleus for the pair production on isotropic photons is given by

$$\frac{d\gamma }{dt}=\alpha r_0^2cZ^2\frac{m_e}{m_A}_2^{\mathrm{}}𝑑\kappa n\left(\frac{\kappa }{2\gamma }\right)\frac{\phi (\kappa )}{\kappa ^2},$$
(1)
where $`\gamma `$ is the Lorenz factor of the particle, $`\kappa =2k\gamma `$ ( $`k`$ is the momentum of the particle in units of $`m_ec`$), $`n(\kappa )`$, $`\alpha `$, $`r_0`$, $`Z`$, and $`m_A`$ are the photon distribution in the momentum space, the fine-structure constant, the classical electron radius, the charge of the particle, and the rest mass of the particle, respectively. For the energy range $`E10^{19.0}`$ eV, $`\phi `$ can be well represented as

$$\phi \kappa \underset{i=0}{\overset{3}{}}d_i\mathrm{ln}^i\kappa ,$$
(2)
$$d_086.07,d_150.96,d_214.45,d_3=8/3.$$
(3)
Since protons lose a large fraction of their energy in the photopion production, its treatment is very important. We use the interaction length and the energy distribution of final protons as a function of initial proton energy which is calculated by simulating the photopion production with the event generator SOPHIA (Mucke et al., 2000). We checked that arrival spectrum at Earth for mono-energetic injection of energy $`E=10^{21.5}`$ eV is consistent with that of Stanev et al.(2000) (Fig. 2), which is also calculated with the event generator SOPHIA, for the case of ($`B`$, $`l_\mathrm{c}`$) = (1nG, 1Mpc) which is the closest parameter set to the one of Stanev et al.(2000).

The EGMF are little known theoretically and observationally. The upper limit for its strength and correlation length $`B_{\mathrm{IGM}}l_\mathrm{c}^{1/2}<1\mathrm{n}\mathrm{G}(1\mathrm{M}\mathrm{p}\mathrm{c})^{1/2}`$, as measured by Faraday rotation of radio signals from distant quasars (Kronberg, 1994), is often used (Stanev et al., 2000; Smialkowski et al., 2002). We adopt in this study not only this widely accepted value, 1nG, but also relatively strong magnetic field (10nG, 100nG), which is better for explaining the extension of the energy spectrum and the observed isotropic arrival distribution. In this case, deflection angles of UHECRs become so large that the no counterparts problem can be simply solved.

We assume that the magnetic field is represented as the Gaussian random field with zero mean and a power-law spectrum. Thus, $`B^2(k)`$ can be written as

$`B^2(k)`$ $`k^{n_H}\mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}2}\pi /l_\mathrm{c}k2\pi /l_{\text{cut}}`$ (5)
$`=0\mathrm{otherwise},`$
where $`l_\mathrm{c}`$ is the correlation length of the EGMF and $`l_{\text{cut}}`$ characterizes the numerical cut-off scale. We use $`n_H=11/3`$ corresponding to the Kolmogorov spectrum. According to our previous study (Ide et al., 2001), the correlation length is chosen to be 1 Mpc, 10 Mpc, and 40 Mpc. 1 Mpc is roughly equal to mean separation of galaxies and widely used as the typical value of the correlation length. As we show later, the small scale clustering observed by the AGASA is difficult to be reproduced for $`l_\mathrm{c}10`$ Mpc (See Figure 3.3 and 3.4). We do not have to consider the correlation length larger than 40 Mpc. Physically one expect $`l_{\text{cut}}l_\mathrm{c}`$, but we set $`l_{\text{cut}}=1/8\times l_\mathrm{c}`$ in order to save the CPU time.

The universe is covered with cubes of side $`l_\mathrm{c}`$. For each of the cubes, Fourier components of the EGMF are dialed on a cubic cell in wave number space, whose side is $`2\pi /l_\mathrm{c}`$, with random phases according to the Kolmogorov spectrum, and then Fourier transformed onto the corresponding cubic cell in real space. We create the EGMF of 20 $`\times `$ 20 $`\times `$ 20 cubes of side $`l_\mathrm{c}`$, and outside it, adopt the periodic boundary condition in order to reduce storage data for magnetic field components. Similar methods for the turbulent magnetic fields have been adopted (Sigl et al., 1999; Lemoine et al., 1999; Isola et al., 2002).

In this study, we assume that the source distribution of UHECRs is proportional to that of the galaxies. We use the realistic data from the ORS (Santiago et al., 1995) galaxy catalog, which is nearly full sky survey, and contains two subcatalogs, one complete to a B magnitude of 14.5 (ORS-m) and the other complete to a B major axis diameter of $`1.^{}9`$ (ORS-d). Another nearly full sky catalog which is often used in this kind of study is the IRAS PSCz Survey (Blanton et al., 2001; Smialkowski et al., 2002). However, it is suspected that some galaxies in the LSC are excluded. Since local overdensity of UHECR sources makes the GZK cutoff less sharp or eliminate it, we use the ORS as a source model rather than the IRAS.

In any statistical acceleration mechanism, particles must be kept confined within the acceleration site. Luminous galaxies are generally larger than dwarf ones, and thus they are expected to be able to accelerate particles to higher energy. Furthermore, it is less well known that luminous galaxies in the LSC distribute outward than faint galaxies, contrary to general clusters of galaxies (Yoshiguchi et al., 2002). Accordingly, we inject UHECRs from galaxies with different limiting magnitudes, and explore the source model which reproduces the current observation. As an example, this selection of luminous galaxies corresponds to restricting the sources of UHECRs to host galaxies of AGNs. It is still unknown how much an ultimate UHECR source contribute to the observed cosmic ray flux. We thus consider the two cases in which all galaxies are the same, and they inject cosmic rays proportional to their absolute luminosity. This allow us to investigate effects of the uncertainty of the cosmic ray luminosity on the results.

Since we consider the galaxies with different absolute magnitudes, we use the magnitude-limited subcatalog, i.e., ORS-m, which is also composed of two subsamples, ESO-m (2437 galaxies) and UGC-m (3279 galaxies). We further restrict ourselves to galaxies within 8000 km s-1 and with $`M_B5`$log$`{}_{10}{}^{}h<0.55`$log$`{}_{10}{}^{}v`$, which corresponds to the apparent magnitude limit of $`m_B=14.5`$, where $`h`$ is the dimensionless Hubble constant $`h=H_0/(100`$ km s-1 Mpc$`{}_{}{}^{1})`$ and $`v`$ is the recession velocity. In the following, we abbreviate $`M_B5`$log$`{}_{10}{}^{}h`$ to $`M_B`$. We exclude very faint nearby galaxies with $`M_B>14.0`$. Galaxies in the zone of avoidance ($`|b|<20^{}`$) are not observed. Our final ORS sample contains 5178 galaxies with 2346 galaxies in the ESO-m and 2832 galaxies in the UGC-m. Figure 1 shows the sky distribution of our ORS galaxies in the equatorial coordinate. The strip between the ESO and UGC region $`(17.5^{}\delta 2.5^{})`$ is covered by the ESGC, which is a part of the diameter-limited subcatalog, ORS-d. The absolute magnitudes are shown in Figure 2.1 as a function of the recession velocity. Throughout the paper, we assume that there are no departures from uniform Hubble expansion for simplicity, and use $`h=0.75`$ in calculating distance between galaxies in the sample and our Galaxy from recession velocity.

In order to calculate the energy spectrum and the distribution of arrival directions of UHECRs realistically, there are two key elements of the galaxy sample to be corrected. First, galaxies in a given magnitude-limited sample are biased tracers of matter distribution because of the flux limit. For each galaxy in the sample, we compensate galaxies which is not included in the sample using the selection function given in Santiago et al.(1996). The positions of compensated galaxies are determined according to Gaussian distribution whose mean is the position of the original galaxy and rms is such that $`\sigma =1/3\times (n\varphi (r))^{1/3}`$, where $`r`$ is the distance from our Galaxy and $`n`$, $`\varphi (r)`$ are number density of galaxies which is also calculated using the selection function, value of the selection function at distance $`r`$ from us respectively. For the reason of the selection effect, the ORS is only sampling the universe out to 8000 km s-1 (see Figure 2.1). We assume that contributions from sources outside it is completely isotropic, and calculate its amount from the number of galaxies inside it. Second, our ORS sample does not include galaxies in the zone of avoidance ($`|b|<20^{}`$) and the ESGC region $`(17.5^{}\delta 2.5^{})`$. In the same way, we assume that the source distribution in this region is homogeneous, and calculate its number density from the number of galaxies in the observed region, which is corrected for the selection effect. In Figure 2.1, we show the sky distribution of the ORS galaxies corrected in the manner explained above.

Finally, we now explain how the angular probability distribution of UHECRs are calculated. At first, protons with a flat energy spectrum are injected isotropically at a given point in each of the 26 energy bins, and then propagated in the EGMF over $`1`$ Gpc for $`15`$ Gyr. Weighted with a factor corresponding to a $`E^2`$ power law spectrum, this provides distribution of energy, deflection angle, and time delay of UHECRs as a function of the distance from the initial point. With this distribution, we can calculate the flux and the arrival directions of UHECRs injected at a single UHECR source. Then, summing contributions from all the sources, we obtain the angular probability distributions of UHECRs. In this paper, we use the distribution of energies and deflection angles integrated over the time delay, assuming that the cosmic ray flux at the earth is stationary. In Figure 2.1, we show number of protons with energies of $`(10^{19.6}10^{20.3})`$ eV integrated over deflection angle, time delay as a function of propagation distance. It is noted that this number is proportional to the cosmic ray flux in this energy range from sources in concentric shells with an infinitely small width, assuming that a source distribution is homogeneous.

### 2.2 Statistical Methods

In this subsection, we explain the two statistical quantities, the harmonics analysis for large scale anisotropy (Hayashida et al., 1999) and the two point correlation function for small scale anisotropy. The harmonic analysis to the right ascension distribution of events is the conventional method to search for global anisotropy of cosmic ray arrival distribution. For a ground-based detector like the AGASA, the almost uniform observation in right ascension is expected. The $`m`$-th harmonic amplitude $`r`$ is determined by fitting the distribution to a sine wave with period $`2\pi /m`$. For a sample of $`n`$ measurements of phase, $`\varphi _1`$, $`\varphi _2`$, $`\mathrm{}`$, $`\varphi _n`$ (0 $`\varphi _i2\pi `$), it is expressed as

$$r=(a^2+b^2)^{1/2}$$
(6)
where, $`a=\frac{2}{n}\mathrm{\Sigma }_{i=1}^n\mathrm{cos}m\varphi _i`$, $`b=\frac{2}{n}\mathrm{\Sigma }_{i=1}^n\mathrm{sin}m\varphi _i`$. We calculate the first and second amplitude from a set of events generated according to predicted probability density distribution of arrival directions of UHECRs.

If events with total number $`n`$ are uniformly distributed in right ascension, the chance probability of observing the amplitude $`r`$ is given by,

$$P=\mathrm{exp}(k),$$
(7)
where

$$k=nr^2/4.$$
(8)
The AGASA 57 events is consistent with isotropic source distribution within 90 $`\%`$ confidence level as shown later. Accordingly, we adopt the amplitude such that $`P=0.1`$ as observational constraint, rather than the most probable amplitude calculated from the observed 57 events.

The two point correlation function $`N(\theta )`$ contains information on the small scale anisotropy. We start from a set of generated events or the actual AGASA data. For each event, we divide the sphere into concentric bins of angular size $`\mathrm{\Delta }\theta `$, and count the number of events falling into each bin. We then divide it by the solid angle of the corresponding bin, that is,

$`N(\theta )={\displaystyle \frac{1}{2\pi |\mathrm{cos}\theta \mathrm{cos}(\theta +\mathrm{\Delta }\theta )|}}{\displaystyle \underset{\theta \varphi \theta +\mathrm{\Delta }\theta }{}}1[\mathrm{sr}^1],`$ (9)
where $`\varphi `$ denotes the separation angle of the two events. $`\mathrm{\Delta }\theta `$ is taken to be $`1^{}`$ in this analysis. The AGASA data shows strong correlation at small angle $`(2^{})`$ with 5 $`\sigma `$ significance of deviation from an isotropic distribution as shown later.

In order to quantify the statistical significance of deviations between models and data we introduce $`\chi _{\theta _{\mathrm{max}}}`$, which is defined by,

$`\chi _{\theta _{\mathrm{max}}}={\displaystyle \frac{1}{\theta _{\mathrm{max}}}}\sqrt{{\displaystyle \underset{\theta _\mathrm{i}=0}{\overset{\theta _{\mathrm{max}}}{}}}{\displaystyle \frac{\left(N(\theta _\mathrm{i})N_{\mathrm{obs}}(\theta _\mathrm{i})\right)^2}{\sigma _\mathrm{i}^2}}},`$ (10)
where $`N_{\mathrm{obs}}(\theta _\mathrm{i})`$ is the value of the two point correlation function obtained by the AGASA observation at angle $`\theta _\mathrm{i}`$, and $`\sigma _\mathrm{i}`$ is statistical error of the $`N(\theta _\mathrm{i})`$ due to the finite number of simulated events, which is set to be the observed one. In this study, we take $`\theta _{\mathrm{max}}`$ to be 10 deg for properly quantifying deviations of models from the sudden increase of correlation at the small angle scale observed by the AGASA.

We compare the model predictions with the existing AGASA data throughout the paper. Therefore, we set the number of simulated events to be that of the AGASA data above $`10^{19.6}`$ eV (57), and restrict the arrival directions of UHECRs in the range of $`10^{}\delta 80^{}`$ when calculating these statistical quantities. There are also the present working or development of large-aperture new detectors, such as the HiRes (Wilkinson et al., 1999) and South and North Auger (Capelle et al., 1998). If they report the arrival directions of observed events, we can make comparisons between the model predictions and their data using the same statistical methods.

## 3 RESULTS

### 3.1 Angular Images of UHECRs

In this subsection, we present the results of the numerical simulations for the angular images of UHECRs in a variety of conditions concerning strength and correlation length of the EGMF and the limiting magnitudes of galaxies. We show the angular images for two energy range $`(10^{19.6}10^{20})`$ eV and above $`10^{20}`$ eV, about at which the GZK cutoff of the energy spectrum is predicted. For comparison, we also show the AGASA events with the angular images. It is noted that the current AGASA data set contains 49 events with energies of $`(10^{19.6}10^{20})`$ eV and 8 events above $`10^{20}`$ eV.

The calculated angular images of UHECRs originating the galaxies more luminous than $`M_{\mathrm{lim}}=14.0`$ (upper), $`18.0`$ (middle), $`20.0`$ (bottom) for number (left) and luminosity (right) weighted sources in the case of $`(B,l_\mathrm{c})=(1,1)`$ are shown in Figure 2 for energy range $`(10^{19.6}10^{20})`$ eV. The AGASA 49 events in this energy range are also shown as circles of radius proportional to their energies. Spatial structure of the LSC is reflected by prominent high intensity region running north and south at $`\alpha 180^{}`$ for $`M_{\mathrm{lim}}=14.0`$ and $`18.0`$. Since the luminous galaxies which distribute outward than the faint ones make a substantial contribution in the case of luminosity weighted sources, the region of high intensity becomes to be less obvious in the case of $`M_{\mathrm{lim}}=14.0`$. In other cases, there is no significant difference simply due to small variance of absolute magnitudes of the individual galaxy. We can see that the angular image becomes to be isotropic and show good correlation with AGASA events as restricting sources to more luminous galaxies. Such correlation of the AGASA events with angular image of UHECRs is also pointed out in the case of galaxies with huge infrared luminosity (Smialkowski et al., 2002).

Even for $`M_{\mathrm{lim}}=20.0`$, there are high intensity regions around $`(\alpha ,\delta )(180^{},10^{})`$ which do not correlate with the AGASA events. However, if we select sources, which contribute to the currently observed cosmic ray flux, from our ORS sample, these high intensity regions can be eliminated. An example in the case of $`M_{\mathrm{lim}}=20.5`$ is shown in Figure 3.4 as we explain later.

Figure 3 shows the predicted angular image of UHECRs for $`(B,l_\mathrm{c})=`$(1,1) (upper), (1,10) (middle), (10,1) (bottom) and $`M_{\mathrm{lim}}=14.0`$ (left), $`20.0`$ (right) with energies of $`(10^{19.6}10^{20})`$ eV. The angular image is distorted for $`(B,l_\mathrm{c})=`$(1,10), (10,1), and dependence on the correlation length is relatively weak because the deflection angle of UHECRs is proportional to $`Bl_\mathrm{c}^{1/2}`$. As we will show later, this distortion of the angular image is crucial for the small scale clustering.

In Figure 4, we show the angular image for $`(B,l_\mathrm{c})=(1,1)`$ (upper), $`(10,1)`$ (middle), $`(100,1)`$ (bottom) and $`M_{\mathrm{lim}}=14.0`$ (left), $`20.0`$ (right) with the energy range $`(10^{20}10^{20.3})`$ eV. (Highest energy of the events observed by the AGASA is $`10^{20.3}`$ eV.) Because protons with energies above $`10^{20}`$ eV must originate within about 50 Mpc from us due to photopion production, only nearby galaxies are able to make a substantial contribution to the angular image so that the spatial structure of the LSC is clearly visible. There is no correlation for $`M_{\mathrm{lim}}=14.0`$ between the AGASA events and the angular image of UHECRs. For $`M_{\mathrm{lim}}=20.0`$, arrival distribution is relatively isotropic, but still has no correlation with the AGASA events. As we conclude later, the AGASA events above $`10^{20}`$ eV might be of different origin, if they are confirmed.

### 3.2 Energy Spectra of UHECRs

We present in this subsection the results of the numerical simulations for the energy spectra of UHECRs. Figure 3.2 shows the energy spectrum in the case of $`B=1`$ (upper), $`10`$ (middle), $`100`$ (bottom) nG and $`M_{\mathrm{lim}}=20.0`$ for number weighted sources. They are normalized arbitrary, but with the same factors for all the values of the correlation length. The substantial deflection in the EGMF results in a significant increase of the path length, then leading to a time delay exceeding the age of the universe or a decrease of energy below $`E=10^{19.6}`$ eV due to pair production for protons injected at sufficiently large distance, as easily seen in Figure 2.1. Thus, the stronger EGMF make the GZK cutoff around $`10^{20}`$ eV less sharp or eliminate it in Figure 3.2.

When the correlation length of the EGMF is small, protons can not be deflected effectively and thus its path length and equivalently energy loss is smaller than that for longer correlation length. This can also be seen in Figure 2.1 and Figure 3.2. Especially, effects of changing the correlation length are large at $`E<10^{20}`$ eV for $`B=10`$ nG and at $`E10^{20}`$ eV for $`B=100`$ nG, where gyroradii of UHECRs are about $`1`$ Mpc and thus deflection is strong enough for the diffusion approximation to become applicable in the case of $`l_\mathrm{c}=10`$ and $`40`$ Mpc.

In order to determine normalization of the energy spectrum and quantify statistical significance of mean deviation between the fitted energy spectrum and the observed one, we introduce $`\chi _{ES}`$ in the similar way as the two point correlation function explained in the previous section. This is defined by

$`\chi _{\mathrm{ES}}={\displaystyle \frac{1}{9}}\sqrt{{\displaystyle \underset{E_\mathrm{i}=10^{19.5}}{\overset{10^{20.3}}{}}}{\displaystyle \frac{\left(dN/dE(E_\mathrm{i})dN_{\mathrm{obs}}/dE(E_\mathrm{i})\right)^2}{\sigma _\mathrm{i}^2}}},`$ (11)
where $`dN/dE_{\mathrm{obs}}(E_\mathrm{i})`$ is the energy spectrum observed by the AGASA at $`E=E_\mathrm{i}`$, and $`\sigma _\mathrm{i}`$ is $`1\sigma `$ error at this energy bin.

The origin of cosmic rays with energies below $`10^{19.5}`$ eV is also one of the major open questions in astro-particle physics. Berezinsky et al.(2002) showed that predicted UHECRs flux from GRBs fall short of the observed flux below $`10^{19.5}`$ eV assuming an injection spectrum $`E^2`$, which is advocated in Vietri (1995) and Waxman (1995), and also used in this study. However, there are another possible extra-galactic and galactic UHECR production sites, such as AGNs and young neutron stars (Blasi et al., 2000). Maximum energy of cosmic ray achieved by them may be lower than $`10^{19.5}`$ eV. In this case, these galactic and/or extra-galactic components may substantially contribute to the cosmic ray flux below $`10^{19.5}`$ eV. Throughout the paper, we assume that these components bridge the gap between the observed flux and the predicted one of Berezinsky et al.(2002), and restrict ourselves to the energy spectrum only above $`10^{19.5}`$ eV.

Figure 3.2 shows the energy spectrum normalized so as to minimize $`\chi _{ES}`$ in the case of $`B=1`$ (upper), $`10`$ (middle), $`100`$ (bottom) nG and $`M_{\mathrm{lim}}=14.0`$ (left), $`20.0`$ (right) for number weighted sources. We also show the significance values $`\chi _{ES}`$ in this figure.

The predicted energy spectrum for $`(B,M_{\mathrm{lim}})=(1,14.0)`$ provide good fit to the AGASA data, including the extension above $`10^{20}`$ eV. However, arrival distribution is significantly anisotropic and inconsistent with the AGASA (see the next subsection). In the case of $`(B,M_{\mathrm{lim}})=(1,20.0)`$, the energy spectrum can be well fitted with the AGASA data below $`E=10^{20}`$ eV, but the GZK cutoff is predicted above $`E=10^{20}`$ eV due to the energy loss by photopion production in contrast to the AGASA data. This is because luminous galaxies near the earth are less than dwarf galaxies (see Figure 2.1). Although the energy spectrum extends beyond $`10^{20}`$ eV for the stronger EGMF, it does not provide good fit with the AGASA data at $`E<10^{20}`$ eV, where the current data have adequate statistics. Consequently, $`\chi _{ES}`$ for the stronger EGMF is larger than that for $`B=1`$ nG as easily seen in Figure 3.2.

In the same way as Figure 3.2, we show the energy spectra for luminosity weighted sources in Figure 3.2. Since contributions from the dwarf galaxies are suppressed in this case, the cutoff of the energy spectra for $`M_{\mathrm{lim}}=14.0`$ is predicted. For $`M_{\mathrm{lim}}=20.0`$, the results are almost the same as number weighted ones.

### 3.3 Statistics on the Arrival Directions of UHECRs

In this subsection, we show the results of statistics on the arrival directions of UHECRs. Figure 3.3 shows the predicted first harmonics for number weighted sources as a function of limiting magnitude. The error bars represent the statistical error due to the finite number of simulated events, whose number is set to be that observed by the AGASA in the energy range of $`(10^{19.6}10^{20.3})`$ eV. Arrival directions of UHECRs are restricted in the range $`10^{}\delta 80^{}`$ in order to compare our results with the AGASA data. The solid lines indicate the harmonic amplitude obtained by the AGASA data. The shaded region is expected from the statistical fluctuation of isotropic source distribution with the chance probability larger than $`10\%`$.

Clearly visible in Figure 3.3 is that the anisotropy is least for $`M_{\mathrm{lim}}20.5`$, which is caused by broadened distribution of luminous galaxies in the LSC. The amplitude at $`M_{\mathrm{lim}}=20.5`$ is consistent with isotropic source distribution within 1 $`\sigma `$ level for all the $`(B,l_\mathrm{c})`$. Compared with the AGASA data, a significant anisotropy can be seen for $`M_{\mathrm{lim}}18.0`$, although the energy spectra in these cases extend above $`10^{20}`$ eV as shown in the previous subsection. Increase of the amplitude around $`M_{\mathrm{lim}}21.0`$ is mainly due to two giant galaxies existing in the almost same direction $`(\alpha 185^{},\delta 14.7^{})`$ in the vicinity of our Galaxy (at $``$ 30 Mpc).

Comparing the results of different strength of the EGMF, we find that the anisotropy in the cases of 10, 100 nG is larger than that of 1 nG below $`M_{\mathrm{lim}}18.0`$ contrary to expectation. We note again that stronger EGMF results in shorter cosmic ray horizon because of time delay and energy loss due to the substantial deflection in the EGMF (See Figure 2.1). That is, spatial structure of nearby galaxies is strongly reflected in the large scale anisotropy for the strong EGMF. Above $`M_{\mathrm{lim}}18.0`$, since luminous galaxies distribute outward, the anisotropy in the cases of the strong EGMF is smaller as expected. Also, increase of correlation length results in more effective deflections, and then leads to smaller anisotropy in the case of 1, 10 nG. However, this dependence is opposite for 100 nG, which is caused by a reason somewhat similar to that of slightly large anisotropy for the strong EGMF below $`M_{\mathrm{lim}}18.0`$ explained above. Gyroradius of proton with energy $`10^{20.0}`$ eV in the EGMF of 100 nG is 1 Mpc, which is smaller than $`l_\mathrm{c}=10,40`$ Mpc. In this case, UHECR deflection is strong enough for the diffusion approximation to become applicable. Thus, the cosmic ray horizon is extremely shorter than that for $`l_\mathrm{c}=1`$ Mpc, so that the anisotropy arise from nearby galaxies in the LSC.

Figure 3.3 shows the predicted second harmonics for number weighted sources in the same manner as the first harmonics. Dependence on the limiting magnitude is weaker than that of the first harmonics. There is no parameter set $`(B,l_\mathrm{c},M_{\mathrm{lim}})`$ which is consistent with isotropic source distribution. We also calculated the first and second harmonics for luminosity weighted sources, and found that the anisotropy is slightly less than that of number weighted sources for $`M_{\mathrm{lim}}18.0`$ because of small contributions from the dwarf galaxies, which concentrate in the LSC. For $`M_{\mathrm{lim}}<18.0`$, the results are almost the same as number weighted sources.

The two point correlation functions are shown for number weighted sources in the case of $`M_{\mathrm{lim}}=14.0`$ (left panels) and $`M_{\mathrm{lim}}=20.5`$ (right panels) in Figure 3.3. The histograms represent the AGASA data, which show the statistically significant correlation at the smallest angle scale $`(1^{})`$. $`\chi _{10}`$ defined in the previous section is also shown.

At larger angles $`(5^{})`$, there are relatively strong correlation for $`M_{\mathrm{lim}}=14.0`$, which is not consistent with the observed isotropic distribution, whereas arrival distributions for $`M_{\mathrm{lim}}=20.5`$ are sufficiently isotropic. Among the cases of $`M_{\mathrm{lim}}=20.5`$, noticeable feature is coexistence of a strong correlation at small angles and sufficiently isotropic distribution at larger scale for $`(B,l_\mathrm{c})=(1,1)`$. This is a consequence of the sharp peak of the UHECR angular image in the case of $`(B,l_\mathrm{c})=(1,1)`$ (see Figure 2). This strong correlation at small angles is reduced or eliminated for longer correlation length or stronger EGMF. Due to small number of observed events, the statistics of the two point correlation function is limited. For this reason, $`\chi _{10}`$ for $`(B,l_\mathrm{c})=(1,1)`$ is equal to that for $`(B,l_\mathrm{c})=(1,10)`$. However, if the future experiments like the Pierre Auger array (Capelle et al., 1998) confirm the small scale clustering of UHECR arrival directions with more event number, the case of $`(B,l_\mathrm{c})=(1,1)`$ will be favored with better significance.

Let us summarize the results of the three quantities (energy spectrum, harmonic amplitude, two point correlation function). First, the energy spectrum for $`B=1`$ nG provide better fit to the AGASA data than that for strong EGMF in terms of $`\chi _{ES}`$. Among this case, the extension of the energy spectrum is also reproduced for $`M_{\mathrm{lim}}14.0`$, although there is significant large scale anisotropy contrary to the AGASA data in this case. In the case of $`M_{\mathrm{lim}}20.0`$, the GZK cutoff is predicted.

Second, the amplitude of the first harmonics is consistent with the AGASA data at $`M_{\mathrm{lim}}20.5`$ irrespective of $`(B,l_\mathrm{c})`$. The amplitude of the second harmonics is not consistent with the AGASA for any $`(B,l_\mathrm{c},M_{\mathrm{lim}})`$.

Finally, the two point correlation function for $`(B,l_\mathrm{c},M_{\mathrm{lim}})=(1,1,20.5)`$ indicate coexistence of a strong correlation at small angles $`(1^{})`$ and sufficiently isotropic distribution at larger scale. In sum, the source model of $`(B,l_\mathrm{c},M_{\mathrm{lim}})=(1,1,20.5)`$ seems to reproduce the observation better than another parameter sets. The amplitude of the second harmonics and the extension of the energy spectrum can not be explained in this source model.

### 3.4 Dependence on the Number of UHECR Sources

Up to present, we have not explained the large scale isotropy quantified by the amplitude of the second harmonics and the small scale clustering except for $`(B,l_\mathrm{c})=(1,1)`$, as well as the extension of the cosmic ray spectrum above $`10^{20}`$ eV. Since we assumed that all the galaxies more luminous than a given limiting magnitude contribute to the cosmic ray flux, the number of sources is maximal and the fluctuations around the assumed non-isotropic distribution is minimal, which make the anisotropy most visible. Furthermore, each source contribute at most one event in this case, it is difficult to obtain clusters which are likely to reflect the point-like sources. Given this situation, we select some galaxies from our ORS sample, and investigate dependence of the results on the number of selected galaxies. It is also expected that we can know the effect of fluctuation of source number within the GZK sphere on the energy spectrum above $`10^{20}`$ eV. As an example, this selection corresponds to restricting the sources of UHECRs to host galaxies of GRBs, which contribute to the currently observed cosmic ray flux.

We assumed that contributions to arrival distribution from sources outside 8000 km s-1 from us, within which the ORS is sampling the universe, is completely isotropic one. However, since we select only some galaxies in the following, point-like nature of sources has to be taken into account. Thus, we distribute galaxies homogeneously in this region, and including them, we select some galaxies from our ORS sample. We restrict ourselves to the limiting magnitude of $`M_{\mathrm{lim}}=20.5`$ and number weighted sources, which produce the most isotropic arrival distribution of UHECRs. It is noted that the number density of galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ in our ORS sample is $`10^{4.3}`$ Mpc-3.

In Figure 3.4, we show the amplitude of the first harmonics as a function of the number fraction (NF) of selected sources to all the ORS galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ for all the $`(B,l_\mathrm{c})`$. For each number fraction, we plot the average over all trial of source selection and realization from the simulated probability distribution with two error bars. The smaller is the statistical error due to the finite number of observed events, while the larger is both the statistical error and the cosmic variance, that is, variation between different selections of sources from our ORS sample. In order to obtain the average and variance, we dial the simulated sets of events 100 times from probability distribution predicted by a specific source distribution, and the sets of sources 30 and 10 times from our ORS sample for the number fraction of $`<0.1`$ and $`>0.1`$ respectively in the case of $`B=1,10`$ nG, but 30 times for all the number fractions in $`B=100`$ nG. We note that the number of sources, which contributes to the cosmic ray flux, for each set of $`(B,l_\mathrm{c})`$ differ from the others, even for the same number fraction, because of the different range of UHECRs (See Figure 2.1). For example, there is no source within 100 Mpc for the number fraction of $`<10^{2.5}`$ (for the number density of $`<10^7`$ Mpc-3) in the case of $`B=100`$ nG.

As a general trend, anisotropy increases with decreasing the number fraction, because the arrival directions of UHECRs clumps to that of sources which is small in number. For large number fraction $`(1)`$, the amplitudes approach to the values in the case that all galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ contribute to the observed cosmic ray flux, as they should. As clearly visible in Figure 2.1, UHECRs can propagate the longest distance for $`(B,l_\mathrm{c})=(1,1)`$, in which the deflection angle, and then the effective path length, is smallest among for another parameters. Accordingly, distant sources make a substantial contribution to the arrival distribution of UHECRs, unlike in the case of stronger EGMF. This is reflected in the decrease of the amplitude around the number fraction of $`10^{1.7}`$ for $`(B,l_\mathrm{c})=(1,1)`$. Except for $`B=100`$ nG, arrival distribution is sufficiently isotropic compared to isotropic source distribution for almost all the number fractions.

The amplitude of the second harmonics is also shown in Figure 3.4. Dependence on the number fraction is roughly same as that of the first harmonics. For $`B=1`$ and $`10`$ nG, the predicted amplitude is sufficiently isotropic for the number fraction of $`(10^{3.0}10^{1.5})`$. On the other hand, it appears to be inconsistent with the observation for any number fraction in the case of $`B=100`$ nG.

Figure 3.4 shows $`\chi _{10}`$ of the two point correlation function as a function of the number fraction. Increase of deviations between models and data with decreasing the number fraction is attributed to the clump of UHECRs at directions of their sources, similarly to the case of the harmonic analysis. Especially, this dependence is quite noticeable for $`(B,l_\mathrm{c})=(1,1)`$, where UHECRs hardly deflect so that correlation at small angle is too strong. At larger number fraction ($`10^{2.0}`$), it can be seen that the model of $`(B,l_\mathrm{c})=(1,1)`$ provide relatively better fit to the observation than another models.

In Figure 3.4, we show the two point correlation function predicted by a specific source scenario for $`(B,l_\mathrm{c})=(1,1)`$, $`(1,10)`$ and $`(10,1)`$ with the number fraction of $`10^{1.7}`$. $`\chi _{10}`$ is also shown. Each source distribution is selected so that it predicts the smallest value of $`\chi _{10}`$ for respective set of $`(B,l_\mathrm{c})`$. Only for $`(B,l_\mathrm{c})=(1,1)`$, another condition is also imposed that the harmonic amplitude predicted by the source model is consistent with the AGASA observation within $`1\sigma `$ level. The correlation of events at small angle scale and sufficient isotropic distribution at larger scale are well reproduced for $`(B,l_\mathrm{c})=(1,1)`$. There is no significant correlation for $`(B,l_\mathrm{c})=(1,10)`$ and $`(10,1)`$, because of the larger deflection angle of UHECRs. In the cases of strong EGMF or longer correlation length, it appears to be difficult to reproduce the strong anisotropy only at the smallest angle scale, even for the case of minimal value of $`\chi _{10}`$.

Considering the calculated results of the first and second harmonics and the small scale anisotropy all together, the model of $`(B,l_\mathrm{c})=(1,1)`$ and the number fraction of $`10^{1.7}`$ (the source number density of $`10^6`$ Mpc-3) seems to reproduce the observations better than another parameter sets. However, due to the small number of observed events, the statistics for the two point correlation function is limited. For this reason, $`\chi _{10}`$ does not differ very much from each other. The future experiments like the Pierre Auger array (Capelle et al., 1998) will decrease statistical uncertainty and provide more strong constraints to the model predictions.

Next we discuss dependence of the energy spectrum on the number fraction. In Figure 3.4, we show the predicted energy spectra in the case of $`(B,l_\mathrm{c})=(1,1)`$ for the number fraction $`10^{0.7}`$, $`10^{1.7}`$, $`10^{2.7}`$ in descending order. They are normalized so as to minimize $`\chi _{ES}`$. The shaded regions represent 1 $`\sigma `$ error due to the source selection from our ORS galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$. As is evident from this figure, decreasing the number fraction increases fluctuation of the predicted energy spectrum, especially above $`10^{20}`$ eV. For the number fraction of $`10^{2.7}`$ (the source number density of $`10^7`$ Mpc-3), since there are no source within the GZK sphere, the GZK cutoff is clearly visible. The predicted energy spectrum is roughly consistent with the AGASA observation within 1 $`\sigma `$ level in the case of the number fraction of $`10^{1.7}`$, where the observed isotropy and clusters are also reproduced well.

However, this result should be interpreted with care. Since the source number density is $`10^6`$ Mpc-3 in the case of the number fraction of $`10^{1.7}`$, mean number of sources is $`0.5`$ within the GZK sphere. Therefore, the AGASA 8 events above $`10^{20}`$ eV, which do not constitute the clustered events with each other, must originate from at most a few sources. We found a source model which explain the extension of the energy spectrum. In this source model, there are two sources within the region covered by the ORS ($`107`$ Mpc), one at $`20`$ Mpc, the other at $`105`$ Mpc. In Figure 3.4, we show the energy spectrum predicted by this source model. The dotted line represent contribution from the nearest source. The predicted energy spectrum provide good fit to the observed one including the extension beyond $`10^{20.0}`$ eV. However, almost all the events above $`10^{20.0}`$ are generated by the nearest source as evident from Figure 3.4. These events are strongly concentrated in a single center, which is inconsistent with the AGASA 8 events. It seems to be difficult to reproduce the AGASA observation above $`10^{19.6}`$ eV (including above $`10^{20}`$ eV) by a single scenario of UHECR origin, even for $`(B,l_\mathrm{c},M_{\mathrm{lim}},NF)=(1,1,20.5,10^{1.7})`$.

Recently, the HiRes Collaboration (Wilkinson et al., 1999) reported the cosmic ray flux from $`2\times 10^{17}`$ eV to over $`10^{20}`$ eV including the GZK cutoff (Abu-Zayyad et al., 2002). In our numerical simulations, the sharp GZK cutoff of the energy spectrum is predicted for the number fraction of $`10^{2.7}`$ as shown above. The small scale correlation of events in this case is too strong for $`(B,l_\mathrm{c})=(1,1)`$ as is evident from Figure 3.4, where the two point correlation function is shown in the same way as Figure 3.4, but with the number fraction of $`10^{2.7}`$. This strong correlation is reduced for $`(B,l_\mathrm{c})=(1,10)`$ and $`(10,1)`$, but still inconsistent with the AGASA data.

Finally we show in Figure 3.4 realizations of 57 arrival directions of UHECRs with energies of $`(10^{19.6}10^{20.3})`$ eV, predicted by the source model of $`(B,l_\mathrm{c},M_{\mathrm{lim}},NF)=(1,1,20.5,10^{1.7})`$ of Figure 3.4. The AGASA 57 events in this energy range are also shown. The arrival directions are restricted in the range $`10^{}\delta 80^{}`$ in order to compare our results with the AGASA data. The arrival distribution of UHECRs seems to be isotropic on a large scale with clusters of events meaning the small scale anisotropy. The amplitudes of the first and second harmonics are consistent with the AGASA, but the extension of the spectrum is not. (Mean number of simulated events above $`10^{20}`$ eV is only 1.) In this example, C6 doublet at $`(\alpha ,\delta )(210^{},35^{})`$ is well reproduced. Also, high intensity regions around $`(\alpha ,\delta )(180^{},10^{})`$ mentioned in Section 3.1. are eliminated under the favor of the source selection.

## 4 DISCUSSION

### 4.1 Constraint on UHECR source model

We presented the results of numerical simulations on the propagation of UHE protons with energies of $`(10^{19.5}10^{22})`$ eV injected by a discrete distribution of galaxies using the ORS data sample, which is properly corrected for the selection effect and absence of galaxies in the zone of avoidance ($`|b|<20^{}`$). We adopted strength of the EGMF to be not only 1 nG but also 10,100 nG, in order to explore the possibility of reproducing the extension of the energy spectrum above $`10^{20}`$ eV. We calculated the three observable quantities, cosmic ray spectrum, harmonic amplitude, and two point correlation function. With these quantities, we made comparisons between the model predictions and the existing AGASA data.

At first, we explored the source model consistent with the current observations, as a function of the limiting magnitudes of galaxies. The energy spectrum for the strong EGMF shows the extension beyond $`10^{20}`$ eV. However, it can not be well fitted with the observed energy spectrum, especially at $`E<10^{20}`$ eV where the current data have adequate statistics. On the other hand, the energy spectrum for $`B=1`$ nG provide good fit with the observation below $`E=10^{20}`$ eV, and statistical significance of deviation $`\chi _{ES}`$ is smaller than that for the strong EGMF. Among this case, the extension of the energy spectrum is also reproduced for $`M_{\mathrm{lim}}14.0`$, although there is significant large scale anisotropy contrary to the AGASA data in this case. In the case of $`M_{\mathrm{lim}}20.0`$, the GZK cutoff is predicted.

We found that the angular image became to be isotropic and showed good correlation with AGASA events as restricting sources to more luminous galaxies. As a consequence, galaxies more luminous than $`20.5`$ mag produce the angular image of UHECR which is most isotropic roughly irrespective of strength and correlation length of the EGMF. However, it is not isotropic enough to be consistent with the observed amplitude of the second harmonics, even for $`M_{\mathrm{lim}}=20.5`$. We further found that the small scale anisotropy is well reproduced in the case of $`(B,l_\mathrm{c})=(1,1)`$ because of the small deflection angle of UHECRs.

Next, in order to obtain sufficiently isotropic arrival distribution of UHECRs and small scale clustering even for the case of other than $`(B,l_\mathrm{c})=(1,1)`$, we randomly selected sources from our ORS sample, and investigate the dependence of the results on their number. We found that the arrival distribution became to be sufficiently isotropic in the case that the number fraction $`(10^{3.0}10^{1.5})`$ of the ORS galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ is selected as UHECR sources for $`B=1`$, $`10`$ nG. Among this case, the small scale clustering can also be reproduced for $`(B,l_\mathrm{c},M_{\mathrm{lim}},NF)=(1,1,20.5,10^{1.7})`$ better than the cases of another parameter sets. We furthermore found that the extension of the energy spectrum is explained in this case within roughly 1 $`\sigma `$ error due to source selection from our ORS sample because of the large fluctuation above $`10^{20}`$ eV.

In sum, we showed that the three observable quantities including the GZK cutoff of the energy spectrum can be reproduced in the case that the number fraction $`10^{1.7}`$ of the ORS galaxies more luminous than $`20.5`$ mag is selected as UHECR sources. In terms of the source number density, this constraint corresponds to $`10^6`$ Mpc-3.

### 4.2 Implication for sources of UHECR above $`10^{20}`$ eV

As mentioned in the previous subsection, we showed that the three observable quantities including the GZK cutoff of the energy spectrum can be reproduced in the case that the number fraction $`10^{1.7}`$ of the ORS galaxies more luminous than $`20.5`$ mag is selected as UHECR sources. However, we should interpret this result with care. Since the AGASA 8 events above $`10^{20}`$ eV do not constitute the clustered events with each other, the number of their sources within the GZK sphere must be much larger than the event number (8). On the other hand, the number of sources within the GZK sphere is at most only 1 in the case of the number fraction of $`10^{1.7}`$ (the source number density of $`10^6`$ Mpc-3). It is suspected that cosmic rays above $`10^{20}`$ eV (hereafter EHECRs; extremely-high energy cosmic rays) are strongly concentrated in a single center, even if the extension of the energy spectrum is explained. Therefore, large fraction of EHECRs observed by the AGASA might originate from sources other than ones of UHECR with energies of $`(10^{19.6}10^{20.0})`$ eV. The number density of EHECR source must be much larger than $`10^6`$ Mpc-3.

We can put forward two candidates as such EHECR production site. It is suggested that iron ions from the surfaces of young strongly magnetized neutron stars may be accelerated to $`10^{20}`$ eV through relativistic MHD winds (Blasi et al., 2000). However, it seems to be difficult to obtain sufficiently isotropic arrival directions of EHECRs from sources localized in the Galactic plane by deflection due to the Galactic magnetic field (O’Neill et al., 2001).

The decay of some supermassive particles, which could be produced from TDs, or be certain MSRPs, are also considered as probable EHECRs origin. Decay rate of the supermassive particles required to explain the observed EHECR flux are estimated as $`10^{35}`$ Mpc-3 yr-1 (Bhattacharjee and Sigl, 2000). Multiplying by the current EHECR observing time ($`10`$ yr), we obtain the number density of $`10^{36}`$ Mpc-3, which is sufficiently large to reproduce not clustered AGASA 8 EHECRs. As for the large scale anisotropy, distribution of TDs in the universe would be homogeneous. MSRPs are expected to be clustered in galactic halos, and thus EHECR flux will be dominated by contribution from our own Galactic Halo. Accordingly, arrival distribution would be isotropic enough to be consistent with the current AGASA observation in both cases. Furthermore, these scenarios are not constrained by contribution to the cosmic ray flux at lower energy which must be smaller than observed one, since hard injection spectrum is generally predicted. For the reasons stated above, the top-down scenarios may be favored to explain large fraction of the 8 EHECRs observed by the AGASA.

There is a problem if all the EHECRs observed by the AGASA originate from sources other than those of UHECRs. If we consider that sources of EHECRs and UHECRs precisely differ from each other, there would be no clustered event set which includes both EHECRs and UHECRs, because clustered events are expected to reflect a single point source. However, 2 EHECRs out of 8 observed by the AGASA constitute the clustered events with UHECRs with energies of $`(10^{19.6}10^{20})`$ eV (C1 and C3 in Hayashida et al.(2000)). We have to interpret that the two clusters happen to be produced by chance. Chance probability of more than two accidental clusterings within $`2.5^{}`$ between 49 UHECRs and 8 EHECRs is quite low $`(4\%)`$. The AGASA observation seems to deviate from our prediction with $`2\sigma `$ confidence level.

However, about 1 EHECR out of 8 may originate in the bottom-up scenarios. Indeed, 57 cosmic rays above $`10^{19.6}`$ eV predicted by the source model of Figure 3.4 include an EHECR. In this case, since the source of this single event will be located within the GZK sphere, it is likely that this event constitutes a clustered event set with UHECRs. Chance probability of more than one accidental clustering rise up to $`28\%`$, which is high enough for us to interpret one of the two doublets to be produced merely by chance. Thus, our scenario is consistent with the current AGASA observation within 1 $`\sigma `$ level.

The event correlation between EHECRs and UHECRs would provide an important test of our conclusion with sufficient amount of data from future experiments. For example, when the event number of UHECRs above $`10^{19.6}`$ increases to 10 times as many as that of the current AGASA observation (57), the number of EHECRs which lie within $`2.5^{}`$ from UHECRs by chance becomes $`18.6\pm 4.3(1\sigma )`$. If future experiments observe event number of such EHECRs within this range, our scenario will be supported. Here we note that the source number within the GZK sphere is at most 1 in the case of NF=$`10^{1.7}`$. Therefore the number of event clustering between an EHECR and an UHECR which can be explained by bottom-up scenarios is at most 1.

Another test of our conclusion is distinguishing the chemical composition of primary cosmic rays which initiate air showers in the earth atmosphere. In the bottom-up scenarios, protons or nuclei are primary particles. On the other hand, photons are considered as primary in the top-down scenarios (Bhattacharjee and Sigl, 2000). Information on the chemical composition is mainly provided by the depth of shower maximum for fluorescence observation of the air shower. There are a number of projects utilizing the fluorescence technique, such as HiRes (Wilkinson et al., 1999), South and North Auger (Capelle et al., 1998), OWL (Cline and Stecker, 2000) and EUSO (Benson and Linsley, 1992). These experiment would provide information on the chemical composition and a test of our conclusion.

### 4.3 Comparison with the HiRes cosmic ray spectrum

The HiRes Collaboration (Wilkinson et al., 1999) reported the cosmic ray flux from $`2\times 10^{17}`$ eV to over $`10^{20}`$ eV including the GZK cutoff (Abu-Zayyad et al., 2002). As discussed above, we showed that the three observable quantities except for EHECR events observed by the AGASA can be reproduced in the case that the number fraction $`10^{1.7}`$ of the ORS galaxies more luminous than $`20.5`$ mag is selected as UHECR sources. In Figure 4.3, we show the predicted energy spectrum fitted to the data measured by Hires-I detector in the case of $`(B,l_\mathrm{c},M_{\mathrm{lim}},NF)=(1,1,20.5,10^{1.7})`$. As easily seen, the fit is good as compared with the middle panel of Figure 3.4. If the flux measured by the HiRes Collaboration is correct and observational features about the arrival distribution are same as the AGASA, our source model can explain both the arrival distribution of UHECRs and the flux at the same time. When the large-area detectors like the Pierre Auger array (Capelle et al., 1998) start operating, the statistical error of the observed energy spectrum will decrease. This allows us to draw a conclusion about the origin of EHECRs or whether they exist or not.

### 4.4 Implication for sources of UHECR below $`10^{20}`$ eV

In this subsection we discuss implications of our results, obtained by comparison of the model predictions with the existing AGASA data, for the origin of UHECRs with energies of $`(10^{19.6}10^{20})`$ eV. The number densities of possible UHECR production sites and ratios to our result ($`10^6`$ Mpc-3) are tabulated in Table 1.

Blast waves in AGN jet could accelerate particle to a few tens of EeV (Halzen and Zas, 1997) and similarly for AGN cores. However, due to the interaction of UHECRs with very high radiation fields in and around the central engine of an AGN, these maxima are unlikely to be achieved. This is consistent with our result, because the number density of AGN are much larger than our results, and thus, arrival distribution may be inconsistent with the AGASA if all the AGNs contribute to the cosmic ray flux above $`10^{19.6}`$ eV. Among AGNs, the most promising acceleration sites for UHECRs are the so-called hot-spots of FR-II radio galaxies (Biermann, 1997; Norman et al., 1995). However, their number density seems to be too small, as compared with our result. It is unlikely that FR-II radio galaxies are responsible for all the AGASA 49 events with energies of $`(10^{19.6}10^{20})`$ eV. Another possible UHECRs production sites are BL Lacertae objects (BL Lacs) (Tinyakov and Tkachev, 2001). Their number are larger than FR-II radio galaxies, but still insufficient to reproduce all the 49 events.

GRB can be shown to in principle accelerate protons up to $`10^{20}`$ eV (Waxman, 2000). The effective density of GRB is determined by the GRB rate $`(2\times 10^{10}`$ h3 Mpc-3 yr-1 (Mao and Mo, 1998)) and the typical time delay $`\tau `$ of UHECRs, which depends on the propagation distance from sources to the earth. Since we are now considering UHECRs with energies of $`(10^{19.6}10^{20})`$ eV, we have to include contribution from sources outside the GZK sphere when calculating the typical time delay. We calculated the mean time delay from our numerical data assuming homogeneous GRBs distribution, and obtained $`\tau 10^6`$ yr. This gives the GRB effective density $`10^4`$ Mpc-3. From our result, if about 1 $`\%`$ of GRBs produce UHECRs, arrival distribution predicted by GRBs would be consistent with the AGASA.

The compact dynamo model has been proposed as a natural mechanism for accelerating cosmic rays in dead quasars (Boldt and Ghosh, 1999; Levinson, 2000; Boldt and Loewenstein, 2000). Only 0.2 $`\%`$ of dead quasars are required to contribute to the comic ray flux above $`10^{19.6}`$ eV.

Large fraction of luminous infrared galaxies (LIRGs) are found to be interacting systems, suggesting that collisions and merging processes are responsible for triggering the huge IR light emission (Sanders et al., 1988). In such systems, favorable environments for accelerating cosmic rays to $`10^{20}`$ eV are provided by amplified magnetic fields on the scale of tens kpc resulting from gravitational compression, as well as high relative velocities of galaxies and/or superwinds from multiple supernovae explosions (Cesarsky, 1992). According to Smialkowski et al.(2002), we take the number density of such colliding galaxies as that of galaxies with huge far infrared luminosity $`L_{\mathrm{FIR}}>10^{11}L_{}`$. The number density is roughly equal to our result. If all the galaxies with $`L_{\mathrm{FIR}}>10^{11}L_{}`$ are responsible for UHECRs, predicted arrival distribution would be consistent with the AGASA. Indeed, Smialkowski et al.(2002) showed that arrival distribution of UHECRs predicted by LIRGs is sufficiently isotropic such that statistical tests they used for large scale anisotropy was not conclusive for distinguishing between isotropic and LIRGs source distribution. Furthermore, they showed that the probability of occurring the clustered events are more than 10 times higher for LIRGs source distribution than that for isotopic distribution. Their results well agree with ours.

### 4.5 Consistency with a statistical analysis of event cluster

Our result is supported by a statistical analysis of clustering of UHECRs developed by Dubovsky et al.(2000), which allow an estimate of the minimum number density of sources. They applied their arguments to the 14 events above $`10^{20}`$ eV observed by several experiments. They obtained the minimum number density of sources $`6\times 10^3`$ Mpc-3. Since we are now considering UHECRs below $`10^{20}`$ eV, we have to apply the arguments given in their paper to the cosmic ray events below $`10^{20}`$ eV.

Let us apply the arguments to the 49 events with energies of $`(10^{19.6}10^{20.0})`$ eV observed by the AGASA. According to them, the number density must be larger than

$`h_{}={\displaystyle \frac{1}{4R^3}}{\displaystyle \frac{N_{\mathrm{tot}}^3}{N_{\mathrm{cl}}^2}},`$ (12)
where $`N_{\mathrm{tot}}`$, $`N_{\mathrm{cl}}`$, $`R`$ are the total number of events, the number of events in clusters, e-folding length of cosmic ray luminosity from a single source respectively. The value of $`R`$ can be obtained as $`700`$ Mpc from Figure 2.1. Taking $`N_{\mathrm{tot}}=49`$ and $`N_{\mathrm{cl}}=9`$, Eq.( 12) gives $`h_{}10^6`$ Mpc-3 which is consistent with our result $`10^6`$ Mpc-3.

### 4.6 Strong EGMF in the LSC

Throughout the paper, we assumed that the statistical properties of the EGMF are uniform. However, simple analytical arguments based on magnetic flux freezing, and large scale structure simulations passively including the magnetic field (Kulsrud et al., 1997) demonstrate that the magnetic field is most likely as structured as are the baryons. The EGMF as strong as $`1\mu `$G in sheets and filaments of large scale galaxy distribution, such as in the LSC, are compatible with existing upper limits on Faraday rotation (Ryu et al., 1998; Blasi et al., 1999). It is suspected that the arrival distribution of UHECRs depends on the fields in the immediate environment of the observer.

However, numerical simulations of UHECR propagation in inhomogeneous EGMF over cosmological distances are highly time-consuming. With the assumption of homogeneous EGMF, we have performed numerical simulation of spherically symmetric propagation. Provided that the EGMF is inhomogeneous, the propagation of UHECRs from a single point source has no longer spherical symmetry. As a result, we have to specify the earth position in the universe. We also have to choose the detector (earth) size so small enough for us to accurately calculate the arrival directions. In this case, the number fraction of UHECRs arriving at the earth to injected ones is extremely small. This requires the number of particle to be propagated several orders of magnitudes higher than that used in this study, which takes enormous CPU time.

Although the strong EGMF in the LSC is compatible with the Faraday ratation measures, the cosmic ray observations do not appear to support such local EGMF. Indeed, we showed in this study that small scale clustering can be well reproduced in the case that UHECRs propagate along nearly straight lines. If the local strong EGMF affects the arrival directions of UHECRs, small scale clustering observed by the AGASA will not be obtained.

This difficulty can be seen from the studies by another authors. There are several groups who study the propagation of UHECR in such strong EGMF in the LSC (Sigl et al., 1999; Lemoine et al., 1999; Isola et al., 2002; Sigl, 2002). They assume simplified source distributions which represent the LSC. The strong EGMF of $`1\mu `$G leads to substantial deflections of UHECRs, which are better for explaining the observed isotropic arrival distribution of UHECRs. However, consistency of small scale anisotropy and also large scale isotropy predicted by their scenarios with the AGASA observation is somewhat worse than that predicted by our scenario. The strong EGMF in the LSC seems to be unfavored by the current cosmic ray observation, although it is not ruled out by the Faraday rotation measurements.

Of course, since the event number of UHECRs is currently very small, we can not absolutely rule out such strong EGMF in the LSC. Fulture experiments would increase the significance level of statement about whether such strong EGMF affects the arrival directions of UHECRs or not.

## 5 CONCLUSION

In this section, we present our conclusions obtained from comprehensive study on the bottom-up scenarios. For the origin of EHECRs, we conclude that large fraction of the AGASA 8 events above $`10^{20}`$ eV might originate in the topdown scenarios, or that the cosmic ray flux measured by the HiRes experiment might be better. We also discussed the origin of UHECRs below $`10^{20.0}`$ eV through comparisons between the number density of astrophysical source candidates and our result ($`10^6`$ Mpc-3). At present, we can not conclude which source candidates are responsible for UHECRs below $`10^{20.0}`$ eV. However, our result for the source number density ($`10^6`$ Mpc-3) will be useful when searching for the ultimate UHECR source.

Throughout the paper, we set constraints to the source models so that their predictions reproduce the current AGASA observation within 1 $`\sigma `$ error of the statistical fluctuation due to finite number of the observed events. However, if we relax the constraints to within 2 $`\sigma `$ error, the energy spectrum (Figure 3.2 and 3.2) and the amplitude of the second harmonics (Figure 3.3) can be reproduced in the case of $`(B,M_{\mathrm{lim}})(1,20.5)`$ without the source selection. In this case, since the source number density is $`5.4\times 10^5`$ Mpc-3, there are $`30`$ sources within the GZK sphere. The AGASA events above $`10^{20}`$ eV, which do not constitute the clustered events with each other, may also be reproduced, contrary to the case of the source number density of $`10^6`$ Mpc-3. Thus, we can not absolutely rule out the bottom-up scenarios as the origin of EHECRs.

The statistics for the quantities which we used is limited by the small number of observed events, which make difficult to draw a decisive conclusion on the origin of UHECRs. However, the present working or development of large-aperture new detectors, such as HiRes (Wilkinson et al., 1999) and South and North Auger (Capelle et al., 1998) will considerably decrease the statistical uncertainties. There are furthermore plans for space based air shower detectors such as OWL (Cline and Stecker, 2000) and EUSO (Benson and Linsley, 1992) which would detect about 1000 events per year above $`10^{20}`$ eV (Bhattacharjee and Sigl, 2000). With the method to constrain the source model and their number density developed in this paper, these experiments will reveal the origin of UHECRs with better significance in the very near future.

We would like to thank Dr. T. Stanev for kindly giving us the data of interaction length and inelasticity of UHE protons due to the photopion production with photons of the cosmic microwave background. This research was supported in part by Giants-in-Aid for Scientific Research provided by the Ministry of Education, Science and Culture of Japan through Research Grant No.S14102004 and No.S14079202.