IPMU09-0111
Preprint typeset using LATEX style emulateapj v. 11/12/01

DOES GALACTIC MAGNETIC FIELD DISTURB THE CORRELATION OF THE HIGHEST ENERGYCOSMIC RAYS WITH THEIR SOURCES?
H AJIME TAKAMI 1 AND K ATSUHIKO S ATO 1,2
hajime.takami@ipmu.jp
IPMU09-0111

ABSTRACTThe propagation trajectories of the highest energy cosmic rays can be deflected by Galactic magnetic field(GMF) and expected correlation between their arrival directions and the positions of their sources can be disturbed. In order to explore whether the possible correlation is disturbed or not, we simulate the arrival distributionof the highest energy protons (HEPs) with energies above 6 × 1019 eV taking 4 different GMF models into accountand investigate the cross-correlation between the protons and sources assumed in the simulation. We show thatthe error of cross-correlation function adopted in this study is sufficiently small by accumulating 200 events. We
also find that the correlation is not disturbed largely in many cases after 200 events accumulation and positive
signals of the correlation are significantly expected at angular scale of 3-5◦ . Only in the cases of the northern sky
with axisymmetric spiral structures of GMF, the cross-correlation functions are consistent with no correlation or
have very low significance even if the correlation is positive at small angular scale. Conversely, no observation of
the correlation within this scale implies no field reversal of GMF outside the solar system. Finally, we draw thepossible source directions of recently published data of the Pierre Auger Observatory (PAO) assuming the composition to be purely protons and irons, and demonstrate that the pure-iron assumption may break the correlationof the PAO data with the large-scale structure of local Universe.
Subject headings: cosmic rays — methods: numerical — ISM: magnetic fields — Galaxy: structure
All these works indicate that the origin of HECRs is extragalactic objects related to the distribution of galaxies. Furthermore,in an optimistic view, we can interpret these results as an evidence that intervening magnetic fields are not as strong as information on the directions of their sources disappears like predictions by Dolag et al. (2005); Takami et al. (2006). (but see Sigl
et al. (2004); Kotera & Lemoine (2008); Das et al. (2008);
Ryu et al. (2009) for strong IGMF). This interpretation implies
that much more accumulation of HECR events can unveil the
positions of their sources as strongly concentrated events withina few degree, as demonstrated by Blasi & de Marco (2004)
without considering any intervening magnetic field. Takami &Sato (2008a) showed that such a correlation is also predicted
for the highest energy protons (HEPs) even considering a structured IGMF model, which reflects matter distribution of local
Universe actually observed. We can regard these works as predictions of the PAO results. However, it is still under debate thatastrophysical objects which correlate with the PAO events arereally HECR sources because these objects have weaker activity than theoretically predicted source candidates (Moskalenko
et al. 2009).
Very recently, the PAO reported that significance of the correlation reported in Abraham et al. (2007, 2008) decreased fortheir new data in the same analysis before, though the correlation is still significant (Hague et al. 2009). Thus, we should reconsider the correlation of HECRs and their sources, and the effect of intervening magnetic fields to possible correlation. The
large deflection angles of HECRs could disturb the correlation.
Possibilities are (i) GMF affects the deflections significantly (ii)the main composition of the HECRs is heavier than protons(iii) the IGMF model used in Takami & Sato (2008a) was optimistic, i.e., the deflections predicted by the IGMF model issmaller than that in the real universe. In this study, we focus
on the possibility (i) assuming that the composition of HECRs

1. INTRODUCTION

The origin of the highest energy cosmic rays (HECRs) is
one of the biggest mysteries in astrophysics. Although theoretical candidates of HECR sources have been proposed (for reviews, see Bhattacharjee & Sigl (2000); Torres & Anchordoqui(2004)), we have not found any source yet. A difficulty to findtheir sources is their deflections by intervening magnetic fields,
i.e., Galactic magnetic field (GMF) and intergalactic magnetic
field (IGMF). Since HECRs are charged hadrons, these magnetic fields deflect the trajectories of HECRs from their sources
to the Earth.
Recent progress to construct HECR observatories with larger
effective area has enabled us to accumulate a number of HECR
events and to draw anisotropic distribution of HECR arrival directions (Takeda et al. 1999; Abraham et al. 2007, 2008).
The anisotropy includes information on the positions of theirsources and intervening magnetic fields from the sources to the
Earth. Also, the Pierre Auger Observatory (PAO) found the correlation between the arrival directions of their 27 events above
5.7 × 1019 eV and extragalactic astrophysical objects in local
Universe within z = 0.018 and angular scale of 3.1◦ (Abraham
et al. 2007, 2008). Several authors have also analyzed the PAO
data since publishing their HECR data. George et al. (2008)
showed a significant correlation of the PAO events with hard
X-ray selected active galactic nuclei (AGNs) by using the 2 dimensional generalization of Kolmogorov-Smilnov test. Moregenerally, the correlation with large-scale structure in the local
Universe had been shown (Ghisellini et al. 2008; Kashti &Waxman 2008; Takami et al. 2009). Especially, Takami etal. (2009) showed that the separation angles between the arrival directions of HECRs and their sources are less than ∼ 15◦
by an analysis independent of specific source candidates, whichbounds the deflection angles of HECRs in the local Universe.

1 Institute for the Physics and Mathematics of the Universe, the University of Tokyo, 5-1-5, Kashiwanoha, Kashiwa, Chiba 277-8568, Japan
2 Department of Physics, School of Science and Engineering, Meisei University, 2-1-1, Hodokubo, Hino-shi, Tokyo 191-8506, Japan

1

2

Takami & Satoare protons and IGMF subdominantly affects the deflections of
HEPs. The other possibilities will be discussed subsequently.
Cosmic rays arriving at the Earth are inevitably affected by
GMF. The importance has driven us to consider the propagationof HECRs in the Galactic space (Stanev 1997; Medina-Tancoet al. 1998; Harari et al. 1999; Kalashev et al. 2001; AlvarezMuniz et al. 2002; Tinyakov & Tkachev 2002; Prouza et al.
2003; Yoshiguchi et al. 2004; Tinyakov & Tkachev 2005;
Kachelriess et al. 2007; Takami & Sato 2008b). These works
have mainly investigated their propagation itself in detail andhave shown that predicted deflections of HEPs are less than
∼ 10◦ , the deflection angles highly depend on their arrival directions, and the apertures of HECR detectors to extragalactic
space are distorted. Some of them guessed HECR sources using information on observed events and their trajectories in the
Galactic space. However, they have not predicted possible signals to find their sources from the arrival distribution of HECRs.
Yoshiguchi et al. (2003) simulated the arrival direction distribution of HEPs using plausible source distribution, taking theirpropagation in a GMF model into account, but they focused on
only the auto-correlation of HEPs.
In this study, we focus on how do we see HEP sources in thearrival distribution of HEPs. We consider the cross-correlation
between the arrival directions of HEPs and the positions of theirsources by simulations taking their propagation in magnetized
Galactic space into account. Since we know the positions of
HEP sources in simulations, we can calculate cross-correlation
functions between them. Cross-correlation function we adopt
in this study can also estimate the angular separation scale ofthe correlation automatically.
This paper is laid out as follows: §2 is devoted to explain
GMF models, a HEP source model, the calculation method of
the arrival distribution of protons, and define cross-correlation
function. In §3, the cross-correlation functions of simulated
HEPs with their sources are calculated for several situations,and the results are interpreted and discussed. Finally, the possibilities of (ii) and (iii) are discussed, and we summarize thisstudy in §4.
2. CALCULATION TOOLS

2.1. GMF Models
The structure of GMF is still controversial as briefly described in our previous paper (Takami & Sato 2008b). Thus,
we also adopt 4 different and extreme models of GMFs fairly
used in that work, which were originally proposed by Stanev(1997). Although we adopted a dipole magnetic field for a z
component of GMF in addition to its spiral component related
to the spiral structure of Milky Way matter in Takami & Sato(2008b), we focus on only the spiral field in this study because
there is no direct evidence for the dipole GMF. We briefly introduce the GMF models below.
The radial and azimuthal components of a spiral magnetic
field in the Galactic plane are given by
Br|| = B(r|| , θ) sin p, Bθ = B(r|| , θ) cos p,(1)

where r|| and θ are the distance from the Galactic center andazimuthal angle around the Galactic center, respectively. Notethat θ is defined as increasing clockwise and θ = 0◦ corresponds
to the Galactic center. p is the pitch angle of the spiral field in
the neighborhood of the solar system, set to p = −10◦ . We assume 2 spiral structures: axisymmetric (AS) and bisymmetric(BS), whose field strengths at a point (r|| , θ) are written as´³ b(r|| ) cos θ − β ln r||r0 ´¯ : BS¯ ³
(2)
B(r|| , θ) =¯
 b(r|| ) ¯cos θ − β ln r|| ¯¯ : AS.
r0
Here β ≡ (tan p)−1 = −5.67 and r0 = 10.55 kpc is the galactocentric distance of the location with maximum field£strength at l ¤=0◦ , which can be expressed as r0 = (R⊕ + d) exp −(π/2) tan p ,where R⊕ = 8.5 kpc is the distance of the solar system from the
Galactic center and d = −0.5 kpc is the distance to the nearest
field reversal from the solar system in the BS model. Negative dmeans that the nearest field reversal takes place inside the solar
system. Note that the AS model does not have any field reversals. b(r|| ) is the radial profile of the strength of the magnetic
field, modelled byR⊕b(r|| ) = B0
,
(3)
r||
where B0 = 4.4µG, which corresponds to 1.5µG in the neighborhood of the solar system. In a region around the Galactic
center (r|| < 4 kpc), the field is highly uncertain and therefore
assumed to be constant and equal to its value at r|| = 4 kpc. The
spiral field is assumed to be zero for r|| > 20 kpc.
For a magnetic field in the Galactic halo, we adopt an exponential decreasing model with two scale heights (AlvarezMuniz et al. 2002),½
exp(−|z|): |z| ≤ 0.5 kpc(4)
B(r|| , φ, |z|) = B(r|| , φ)−|z|3exp( 4 − 8 ) : |z| > 0.5 kpc
where the factor exp(−3/8) lets the field continuous on z. The
parities of the spiral fields½ are represented as
B(r|| , φ, z): S type parity
B(r|| , φ, −z) =(5)−B(r|| , φ, z) : A type parity.
The combination between the 2 spiral structures and 2 parities
leads to 4 different GMF models. We adopt these 4 GMF models in this study.
2.2. Source model of HEPs
We consider HEP sources distributed to be compatible withthe large-scale structure (LSS) of local Universe actually observed. The adopted GMF models reflect the structure of our
Galaxy to some extent and the deflection patterns of HEPs are
expected to be dependent on the arrival directions of them. Furthermore, the PAO results imply HEP sources distributed to becomparable with the LSS of local Universe (Abraham et al.
2007, 2008; Ghisellini et al. 2008; Kashti & Waxman 2008;
Takami et al. 2009). Thus, such a source distribution is essential. A model used in this study is originally based on Takami
et al. (2006). Here, the model is briefly explained.
In the model, we regard subsets of galaxies as HEP sources.
We adopt Infrared Astronomical Satellite Point Source Redshift
Survey (IRAS PSCz) catalog of galaxies (Saunders et al. 2000)
to construct the local distribution of galaxies. The IRAS PSCz
catalog matches our purpose well, since it is a flux-limited, uniform catalog and has very large sky coverage (∼ 84% on thewhole sky). However, this catalog includes 2 selection effects:
one is unobserved galaxies with fluxes below the flux-limit of
the IRAS satellite. The other is ∼ 16 % unobserved area whichcorresponds to mainly the Galactic plane. We correct theseselection effects by using a luminosity function of the IRAS
galaxies (Takeuchi et al. 2003). For the former effect, we
add virtual galaxies with luminosities less than the observationally limited flux in the vicinity of original IRAS galaxies. On
the other hand, we simply assume that galaxies are distributed

Does GMF Disturb the Correlation of HECRs with their Sources?

3

3000
BS-S
BS-A
AS-S
AS-A

2500

2000

z [kpc]

1500

1000

500

0

-500

-1000
0

5000

10000
x [kpc]

5000

15000

20000

20000
BS-S
BS-A
AS-S
AS-A

4000

BS-S
BS-A
AS-S
AS-A

15000

3000
10000
2000
5000
z [kpc]z [kpc]

1000
0
-1000

0

-5000

-2000
-10000-3000-15000

-4000
-5000
0

5000

10000

15000
x [kpc]

20000

25000

30000

-20000-20000-15000-10000-5000

0x [kpc]

5000

10000

15000

20000

F IG . 1.— Examples of trajectories of protons with the energy of 1019.8 = 6.3 × 1019 eV projected into a plane perpendicular to the Galactic disk and including
the Galactic center and the Earth. Their arrival directions are (`, b) = (180.0◦ , 0.3◦ ) (upper panel), and (180.0◦ , −0.3◦ ) (lower left panel), and (0.0◦ , −44.8◦ ) (lower
right panel), respectively. There are 4 trajectories of the protons in each panels, which correspond to trajectories in the BS-S (red), BS-A (green), AS-S (blue), and
AS-A (magenta) GMF models.

isotropically for the latter one. Then, we randomly select galaxies according to a given number density of HEP sources, ns , tomake a subset of galaxies and regard it as a HEP source distribution. Such a structured source distribution is adopted up to
200 Mpc from our Galaxy because of the completeness of the
IRAS PSCz catalog. Outside 200 Mpc, the distribution of HEP
sources is assumed to be isotropic. We focus on HEPs above
6 × 1019 eV in this study, 90% of which can reach the Earthfrom their sources within 200 Mpc . Thus, the isotropically
distributed sources does not strongly affect our results.
For the physical properties of HEP sources, we assume that
all the sources emit HEPs persistently with the same power witha power-law spectrum of ∝ E −2.6 whose index can well reproduce the observed spectrum of HECRs above 1019 eV. We consider 2 number densities of HEP sources, ns ∼ 10−4 and 10−5
Mpc−3 . The former can best reproduce the anisotropy of thearrival distribution of HECRs observed by the PAO (Takami &Sato 2009; Cuoco et al. 2008). The latter is the number density
to reproduce the AGASA anisotropy well (Takami et al. 2006;
Blasi & de Marco 2004; Kachelriess & Semikoz 2005) andalso is allowed by the PAO if plausible magnetic fields in the
Universe are taken into account (Takami & Sato 2009).
2.3. Simulation of Arriving HEPs
In order to simulate the arrival distribution of protons fora given source distribution, this study adopts a calculation
method developed by Takami et al. (2006). Generally, it costs
much time to calculate arriving cosmic rays in a magnetized
Universe because whether a cosmic ray injected from a source
to a direction can reach the Earth is not known until its trajectory is calculated. Our method solved this problem and enables

us to calculate the arrival distribution of protons taking theirpropagation in a magnetized space into account in reasonable
CPU time.
In this method, we focus on an inverse process of propagation of protons. The trajectories of protons are calculated by the
backtracking method. The trajectory of an oppositely-charged
proton ejected from the Earth can regarded as that of a proton coming from extragalactic space. Many oppositely-charged
protons (1,250,000 particles per logarithmic bin of energy inthis study) are ejected from the Earth isotropically, and calculate their trajectories in the Galactic space. Once the particles
escape from the Galactic space, they fly away straightforwardly
with energy-gain until the comoving distances of the positionsof the particles exceed 1,500 Mpc or the energy of the particles
exceeds 1022 eV, which corresponds to the maximum injectionenergy of protons at sources, or the propagation time exceeds
the age of the Universe. The boundary between the Galactic and extragalactic space is assumed to be 40 kpc from the
Galactic center. We consider Bethe-Heitler pair creation, photopion production with cosmic microwave background (CMB)
photons and adiabatic energy-loss due to the cosmic expansion
as the energy-loss processes of HEPs in intergalactic space. We
adopt Ωm = 0.3, Ωλ = 0.7, H0 = 71 km s−1 Mpc−1 for the cosmological parameters. For photopion production, we adopt inelasticity and mean-free path calculated by an event generator
SOPHIA (Mucke et al. 2000). The energy-loss is neglected in
the Galactic space because all the energy-loss lengths of protons are much larger than the scale of the Galactic space.
When the trajectory of the ith oppositely-charged proton
passes over sources, it obtains a weight for a positive arriving

4Takami & Sato

F IG . 2.— Arrival directions of HEPs with energies above 1019 eV emitted from sources located at 50 Mpc from the Earth (black circles) in Galactic coordinates.
The energies of HEPs are represented in color. The 4 different GMF models are taken into account.

proton as,
Pselec (E, i) ∝

X
j

j

1
(1 + zi, j )di, j2

dN/dEg (di, j , Eg ) dEg,
E −1.0dE(6)where j labels sources on the trajectory and a factor sums up
for j. zi, j and di, j are the redshift and comoving distance ofthe jth source, respectively. Eg j is the energy of the particle
when it reaches the jth source. dN/dEg (di, j , Eg ) ∝ Eg −2.6 isthe spectrum of proton ejection of oppositely-charged protons
at the jth source. dEg /dE can be calculated by the propagation data. E −1.0 is the appearance of isotropic ejection with
dN/dlog10 E ∝ const. from 1019 to 1021 eV. The weight corresponds to a relative probability in order that the ith particle isa really arriving proton. Thus, we randomly select the required
number of trajectories from all the trajectories according toPselec (E, i), and then arriving protons are obtained. The energy
and arrival directions of the protons are the energy and the injection direction of corresponding oppositely-charged protons.
We take the angular accuracy of HECR observatories to determine the arrival directions of HECRs. We assume that the arrival directions have a fluctuation described by a 2 dimensional
Gaussian distribution with σ = 1◦ , which corresponds to the angular accuracy.
2.4. Statistical Methods
We will discuss the correlation between the arrival directionsof HEPs and their sources in the next section. Cross-correlation
function for this purpose is defined as (Takami et al. 2009;
Blake et al. 2006),
ES(θ) − ES0 (θ) − E 0 S(θ) + E 0 S0 (θ),(7)
wes (θ) =E 0 S0 (θ)
where the characters of E and S represent HEP events and theirsources about which the cross-correlation is calculated, respectively. ES(θ) is the number of pairs between E and G with theseparation angles from θ to θ + ∆θ divided by Ne Ns for normalization, where Ne and Ns are the total numbers of the events andsources considered, respectively. ∆θ is set to be 1◦ . E 0 and
G0 represent HEP events and sources randomly distributed following the apertures of HECR observatories and observationsof HECR sources, respectively. EG0 (θ), E 0 G(θ), and E 0 G0 (θ)
have similar meanings to EG(θ). E 0 and G0 correct their nonuniform apertures. Note that the apertures of HEP sources is
uniform in this study (see §2.2). By definition, w(θ) > 0 means
positive correlation.
The aperture of a ground array is simply described as (Sommers 2001)
ω(δ) ∝ cos(a0 ) cos(δ) sin(αm ) + αm sin(a0 ) sin(δ),(8)
where αm is given by
if ξ > 1 0
πif ξ < −1αm =(9) cos−1 (ξ) otherwiseandcos(θ) − sin(a0 ) sin(δ).
(10)
cos(a0 ) cos(δ)
Here, a0 and θ are the terrestrial latitude of a ground array andthe zenith angle for an experimental cut. We simulate the apertures of 2 HECR observatories; the PAO and Telescope Array(TA). The PAO observes the southern sky and the TA observes
the northern sky. As we will see in §3.1, the deflection patternsof HEPs in the southern sky are quite different from those in
the northern sky. The different signals of correlation between
HECRs and their sources are expected. Thus, we consider theapertures of these observatories as representatives of HECR observatories in the southern and northern sky. The parameters are
a0 = −35.2◦ and θ = 60◦ for the PAO (Abraham et al. 2008),
and a0 = 39.3◦ and θ = 45◦ for the TA (Nonaka et al. 2009).
ξ≡

3. RESULTS

3.1. Deflections of protons by GMF
First of all, we briefly view the trajectories and deflections of
HEPs by GMF. The 4 GMF models predict different trajectories

Does GMF Disturb the Correlation of HECRs with their Sources?

5

F IG . 3.— Deflection angles of arriving protons with energies of 1019.8 eV in color. The arrival directions are shown in equatorial coordinates. The difference
between the AS GMF models and BS GMF models is obvious in the northern hemisphere.

of HEPs, since they are different in field reversals and parities
above and below the Galactic plane.
Fig. 1 shows the trajectories of protons with the energy of1019.8 eV and the same arrival directions in 4 GMF models.
Practically, we eject oppositely-charged protons from the Earthand plot their trajectories. In order to see the effects of different
GMF models, we select 3 examples of the arrival directions ofprotons. In the upper panel, we draw the trajectories of protons
with the arrival direction of (`, b) ∼ (180.0◦ , 0.3◦ ). The trajectory of a proton propagating in the BS-S (AS-S) GMF modelis the same as that in the BS-A (AS-A) GMF model, since aspiral field in the BS-S (AS-S) model is, by definition, the sameas that in the BS-A (AS-A) model above the Galactic plane.
The trajectories of protons are the same both in the S models
and in the A models as long as they propagate above the Galactic plane. The trajectories of the protons in the BS models are
separated from those in the AS models just outside 15 kpc because of a field reversal of GMF. Once a proton passes over a
field reversal, its deflection direction is changed into the opposite direction and the total deflection is weaken. This mechanism works only in the BS models. Thus, generally, the BSmodels predict smaller deflections of HECRs than the AS models. Once oppositely-charged protons propagate away from the
Galactic plane, their deflections become smaller because a magnetic field becomes exponentially smaller. Therefore, the spiralstructure of GMF just in the vicinity of the solar system mainly
affects the deflections of HEPs.
We can see the effect of the parity of GMF in the other panels. The trajectories of protons with the arrival directions of(`, b) = (180.0◦ , −0.3◦ ) are described in the lower left panel. The
trajectories of oppositely-charged protons are deflected toward
the upside of the Galactic plane near the solar system in the caseof GMF models with the S-type parity. On the other hand, theirtrajectories are deflected toward the downside of the Galactic
plane for the A-type parity. Note that the magnetic structure of
the BS-S (BS-A) model is the same as that of the AS-S (AS-A)
model until the first field reversals inside and outside the so-lar system. The trajectories of the oppositely-charged protons
are also separated at ∼ 15 kpc due to a field reversal of GMF.
The lower right panel, in which the arrival directions of HEPs
considered are (`, b) = (0.0◦ , −44.8◦ ), is another example that
HEP trajectories with the same arrival direction are separated
by different GMF models.
The deflection directions of HECRs highly depend on thepositions of their sources, their energies, and of course GMF
structure. Fig. 2 demonstrates the arrival directions of protons
above 1019 eV emitted from 8 extragalactic sources located at
artificial positions shown by black circles to see the directions
and strengths of their deflections. The energies of arriving protons are represented in color and the distance of sources is 50
Mpc.
We can see the dependency of the deflection angles on thepositions of HEPs. For instance, protons from a source with(`, b) ∼ (−60◦ , 60◦ ) are little deflected even for ∼ 1019 eV in
the cases of the AS GMF models. Even for the BS models, theprotons with energies of ∼ 1019 eV are deflected by ∼ 20◦ . On
the other hand, the deflections of protons from a source with(`, b) ∼ (120◦ , 60◦ ) are much larger for their energies of ∼ 1019
eV.
The arrival directions of protons from a source are arranged
reflecting their energies and the GMF structure of the source
direction (Yoshiguchi et al. 2003). The A-type parity leads
to a symmetric structure of deflection pattern of HEPs. This
arrangement structure of the arrival directions of protons canbe a hint to understand the GMF structure of the source directions. When we focus on protons with lower energies, someevents are highly deflected and then their arrival directions
do not reflects the positions of their sources. Furthermore, alarge fraction of such protons is contribution from cosmologically distant sources, and they work as background. So, thelower energy component of the arrangement is hidden. On theother hand, since the highest energy component of protons (redpoints) reaches the Earth only from nearby Universe because ofGreisen-Zatsepin-Kuz’min (GZK) effect (Greisen 1966; Zat-

6

Takami & Sato

Cross-correlation function:wes(θ)

E > 60 EeV, d < 75 Mpc
12120

No GMF
BS-S model
AS-S model

Cross-correlation function:wes(θ)

ns = 10-4 Mpc-3, Np = 27

14

10
8
6
4
2
0

E > 60 EeV, d < 75 Mpc

No GMF
BS-S model
AS-S model

515

ns = 10-5 Mpc-3, Np = 27

100
80
60
40
20
0

-2
0

5

10
Separation Angle:θ [deg]

15

20

0

10
Separation Angle:θ [deg]20

F IG . 4.— Cross-correlation functions between the positions of HEP sources within 75 Mpc and the arrival directions of simulated HEPs above 6 × 1019 eV in thecases of ns ∼ 10−4 (left) and 10−5 Mpc−3 (right). The number of the protons is set to be 27 events and the aperture of the PAO is taken into account. These functions
are calculated under no GMF (green), the BS-S model (red), and the AS-S model (blue).

sepin & Kuz’min 1966), their arrangement is expected to be
more clear, although the deflection angles are small. Such acomponent could be used to investigate GMF structures due totheir small (but not too large) deflections. Also, such small deflections are useful to investigate their sources. Thus, we focus
on the highest energy component of cosmic rays.
Finally, we see the distribution of the deflection angles of
HEPs. in the terrestrial sky. Fig. 3 shows the deflection angles of protons arriving at the Earth with energies of 1019.8 eV
in equatorial coordinates. We can find that GMFs with the AS
structure predict larger deflection than those with the BS structure, as we discussed on Fig. 1. The difference between the AS
models and BS models is noticeable in the northern sky. Thus,
we can predict that the correlation of HEPs and their sources
is quite different between the AS and BS GMF models in thenorthern sky. Note that the large deflection drawn in the lower
right parts of all panels is generated by a strong magnetic fieldnear the Galactic center.
3.2. Cross-correlation with HEP sources
Firstly, we see the effect of GMF to the cross-correlation.
Fig. 4 shows the cross-correlation functions between the positions of sources within 75 Mpc and the arrival directions ofsimulated protons with energies above 6 × 1019 eV in the casesof ns ∼ 10−4 (left) and 10−5 Mpc−3 (right). The number of protons is set to that of the published PAO events (27 events).
Also, the anisotropic aperture of the PAO is taken into account.
These cross-correlation functions are calculated for unmagnetized Galactic space (green), the BS-S GMF model (red), and
the AS-S GMF model (blue). These cross-correlation functions
are calculated as follows; We simulate an arrival distributionof a required number of HEPs for a given source distribution.
The arrival distribution and the source distribution enable us tocalculate cross-correlation function. We realize 100 source distributions with the same ns following the method of §2.2 andcalculate cross-correlation function in the same way for each ofthem. Finally, we calculate the average and standard deviationof the 100 cross-correlation functions and plot them as points
and error bars, respectively. Therefore, the error bars include
not only errors due to the finite number of simulated events but
also errors originating from the uncertainty of HEP source positions, i.e., the uncertainty of galaxy sampling in our HEP source
model. We can find that the correlation signals at small angu-

lar scale are suppressed by GMF. Note that a finite accuracy
to determine the arrival directions of HEPs considered leads topositive values of the cross-correlation functions for no GMF at
the second and third bins.
When we consider a smaller number density of HEP sources,
the absolute values of the cross-correlation functions becomelarger. However, the error bars of the cross-correlation functions also become larger because the uncertainty of source positions is larger. Thus, statistical features are little changed byns .
Next, we simulate a current status of the PAO for 4 GMF
models. Fig. 5 shows the cross-correlation functions of HEP
sources within 75 Mpc and protons above 6 × 1019 eV simulated for the BS-S (red), BS-A (green), AS-S (blue), and AS-A
(magenta). The number of protons is assumed to be 27. We
take the aperture of the PAO into account in the left panel. For
comparison, the aperture of the TA is considered in the rightpanel.
In the left panel, all the cross-correlation functions are notconsistent with a lack of the correlation at a few degree withinstandard deviations, though at the first and second bins these areconsistent with zero because large error bars. Thus, the correlation of HEPs and their sources is expected with a few degree
separation even at present. It is remarkable that the angular separation of the correlation is comparable with the the correlation
scale by the PAO (3.1◦ ).
We can find the difference of cross-correlation signals between GMF models in the northern sky (right panel). The predictions of the AS models are consistent with no correlation,
while the BS models predict positive correlation at a few degree
like the case of the southern sky. This difference is an appearance of the difference of the deflection angles of HEPs between
the AS and BS GMF models, as we have seen in Fig. 3. We
can also see that the 2 BS models (or the 2 AS models) predict
quite similar cross-correlation functions. This means that thedifference between the S-type parity and A-type parity does notappear on the cross-correlation function. Although the difference between the parities is the deflection directions of HEPs
as we saw in Fig. 2, the cross-correlation function does not include information on the deflection directions. These features
will be more clear by accumulating more events as discussed
subsequently.

Does GMF Disturb the Correlation of HECRs with their Sources?
40

ns = 10-5 Mpc-3, Np = 27
E > 60 EeV, d < 75 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

30

Cross-correlation function:wes(θ)

Cross-correlation function:wes(θ)

40

20

10

0

-10

7

ns = 10-5 Mpc-3, Np = 27
E > 60 EeV, d < 75 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

30

20

10

0

-10
0

5

10
Separation Angle:θ [deg]

15

20

0

10
Separation Angle:θ [deg]20

F IG . 5.— Cross-correlation functions of 27 protons for the 4 GMF models with energies above 6 × 1019 eV and their sources within 75 Mpc. ns = 10−5 Mpc−3 isassumed. The anisotropic apertures of the PAO (left) and TA (right) are taken into account.

35

ns = 10-5 Mpc-3, E > 60 EeV
d < 75 Mpc

Cross-correlation function:wes(θ)

30
25

27 events
100 events
200 events
500 events
1000 events

20
15
10
5
0
-5
-10
0

5

10
Separation Angle:θ [deg]1520

F IG . 6.— Cross-correlation functions between HEPs simulated in the BS-S
GMF model with energies above 6 × 1019 eV and their sources within 75 Mpc.
ns = 10−5 Mpc−3 is assumed. The numbers of simulated protons are 27 (red),
100 (green), 200 (blue), 500 (magenta), and 1000 (black), respectively. The
decrease of error bars is saturated at ∼ 200 protons.

The error bars of the cross-correlation functions originate
from not only the finite number of events but also the uncertainty of positions of HECR sources (i.e., the uncertainty ofsampling galaxies in our source model). We can decrease theformer errors by increasing the number of observed events,
whereas the latter type of errors cannot be reduced because we
do not know any source at present. For an ideal case where
the number of detected events is infinite, we can consider only
the latter errors. Practically, it is meaningful to estimate thenumber of events for saturation of the errors because we candiscuss the correlation of HEPs and their sources based on almost only the uncertainty of the precise positions of the sources
after the accumulation of such number of events. Fig. 6 shows
the cross-correlation functions between HEP sources within 75
Mpc and arriving protons simulated for the BS-S GMF modelon the assumption of the PAO aperture. The numbers of protons considered are 27 (red), 100 (green), 200 (blue), 500 (magenta), and 1000 (black). When the number of events is small,errors due to the finite number of events occupy a significant
fraction of the total errors. As increasing the number of events,
the errors decreases and errors due to the sampling of galaxies
are dominated. An estimation by eye implies that the total errors are sufficiently saturated by ∼ 200 events detection. Thisis unchanged for other nS , GMF models, and the northern sky.
Thus, we consider 200 events for discussions below. If we ob-serve the positive value of cross-correlation function after 200
event accumulation, it means that positive correlation between
HEPs and their sources is expected even if the uncertainty of
HEP source positions is taken into account.
Fig. 7 shows the cross-correlation functions of 200 protons
with energies above 6 × 1019 eV and their sources for the BS-S
(red), BS-A (green), AS-S (blue), and AS-A (magenta). The
assumed number density of the sources is ∼ 10−5 Mpc−3 . The
apertures of the PAO (left) and TA (right) are taken into account. We can see the results of Fig. 5 more clearly because of
the smaller errors.
In the left panel, all the 4 GMF models predict the positive
correlation within standard deviations at small angular scale.
Thus, the correlation between HEPs and their sources is predicted despite the uncertainty of the source positions in thesouthern sky. Separation angles of which significance is thehighest is 3◦ and 4◦ for the BS models and the AS models, respectively. Separation angles at which positive correlation ispredicted are smaller in the BS models than in the AS models,reflecting the difference of the deflection angles of HEPs. They
are also up to ∼ 6◦ .
On the other hand, the predictions of the BS models are different from those of the AS models in the northern sky. In
the right panel, the cross-correlation functions calculated for
the AS models are consistent with no correlation within ∼ 3◦within standard deviations, whereas the cross-correlation functions are significantly positive for the BS models at the sameangular scale . At larger angular scale, the cross-correlation
functions are slightly positive for the AS models, but the significance is quite low. Thus, the predictions of the AS GMF models are almost consistent with no correlation. Since the difference between the BS and AS models is based on the existenceof field reversals, observations of very weak or no correlation
at the small angular scale imply no field reversals outside thesolar system.
We can see the difference between the BS and AS GMF
models more obviously by demonstrating the cross-correlation
functions calculated from a given source distribution with ∼10−5 Mpc−3 , shown in Fig. 8. The apertures of the PAO (left)
and TA (right) are considered. We realize the arrival distribution of 200 HEPs 100 times from the source distribution, calculate a cross-correlation function for each of them, and estimate
the average and standard deviation of the 100 cross-correlation
functions. Thus, the error bars originate from only the finite

8

Takami & Sato25
ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 75 Mpc

20

BS-S model
BS-A model
AS-S model
AS-A model

Cross-correlation function:wes(θ)

Cross-correlation function:wes(θ)

25

15

10

5

0

-5

20

ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 75 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

15

10

5

0

-5
0

5

10
Separation Angle:θ [deg]

15

20

0

10
Separation Angle:θ [deg]20

F IG . 7.— Cross-correlation functions between 200 protons simulated for the 4 GMF models with energies above 6 × 1019 eV and the positions of their sources
within 75 Mpc in the case of ∼ 10−5 Mpc−3 . The apertures of the PAO (left) and TA (right) are taken into account.

ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 75 Mpc

15

20
BS-S model
BS-A model
AS-S model
AS-A model

Cross-correlation function:wes(θ)

Cross-correlation function:wes(θ)

20

10

5

0

-5

ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 75 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

15

10

5

0

-5
0

5

10
Separation Angle:θ [deg]

15

20

0

10
Separation Angle:θ [deg]20

F IG . 8.— Cross-correlation functions between 200 protons simulated from a source distribution with ns ∼ 10−5 Mpc−3 with energies above 6 × 1019 eV and thepositions of their sources within 75 Mpc from the Milky Way The apertures of the PAO (left) and TA (right) are taken into account.

5
ns = 10-4 Mpc-3, Np = 200
E > 60 EeV, d < 75 Mpc

4BS-S model
BS-A model
AS-S model
AS-A model

Cross-correlation function:wes(θ)

Cross-correlation function:wes(θ)

5

3

2

1

0

-1

4

ns = 10-4 Mpc-3, Np = 200
E > 60 EeV, d < 75 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

3

2

1

0

-1
0

5

10
Separation Angle:θ [deg]F IG . 9.— Same as Fig. 7, but for ns ∼ 10−4 Mpc−3 .

15

20

0

10
Separation Angle:θ [deg]20

Does GMF Disturb the Correlation of HECRs with their Sources?

3

3

Cross-correlation function:wes(θ)

BS-S model
BS-A model
AS-S model
AS-A model

2.5
Cross-correlation function:wes(θ)

ns = 10-4 Mpc-3, Np = 200
E > 60 EeV, d < 100 Mpc

2.5
2
1.5
1
0.5
0
-0.5

ns = 10-4 Mpc-3, Np = 200
E > 60 EeV, d < 100 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

2
1.5
1
0.5
0
-0.5

-1

-1
0

5

10
Separation Angle:θ [deg]

15

20

0

15

10
Separation Angle:θ [deg]20

15
BS-S model
BS-A model
AS-S model
AS-A model

Cross-correlation function:wes(θ)

ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 100 Mpc
Cross-correlation function:wes(θ)

9

10

5

0

-5

ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 100 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

10

5

0

-5
0

5

10
Separation Angle:θ [deg]

15

20

0

10
Separation Angle:θ [deg]20

F IG . 10.— Cross-correlation functions of 200 protons simulated for the 4 GMF models with energies above ∼ 6 × 1019 eV and the positions of their sources
within 100 Mpc. ns ∼ 10−4 Mpc−3 (upper panels) and ns ∼ 10−5 Mpc−3 (lower panels) are assumed. The anisotropic apertures of the PAO (left panels) and TA(right panels) are taken into account.

10

Takami & Satons = 10-4 Mpc-3, Np = 200
E > 60 EeV, d < 50 Mpc

15
BS-S model
BS-A model
AS-S model
AS-A model

Cross-correlation function:wes(θ)

Cross-correlation function:wes(θ)

15

10

5

0

-5

BS-S model
BS-A model
AS-S model
AS-A model

5

15

10

5

0

-5

80

5

10
Separation Angle:θ [deg]ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 50 Mpc

60

15

20

0
80

BS-S model
BS-A model
AS-S model
AS-A model

Cross-correlation function:wes(θ)

0Cross-correlation function:wes(θ)

ns = 10-4 Mpc-3, Np = 200
E > 60 EeV, d < 50 Mpc

40

20

0

-20

10
Separation Angle:θ [deg]ns = 10-5 Mpc-3, Np = 200
E > 60 EeV, d < 50 Mpc

BS-S model
BS-A model
AS-S model
AS-A model

5

15

60

20

40

20

0

-20
0

5

10
Separation Angle:θ [deg]

15

20

0

10
Separation Angle:θ [deg]20

F IG . 11.— Same as Fig. 10, but for HEP sources within 50 Mpc.

number of events. While all the 4 GMF models show similarcross-correlation functions with significantly positive signals at
small angular scale in the southern sky, cross-correlation functions are clearly different between the AS and BS GMF models in the northern sky. In this example, the BS models predict
strongly positive correlation, but the cross-correlation functions
derived from the AS models are consistent with zero at small
angular scale.
When we consider larger number density of HEP sources,
the values of the cross-correlation function become smaller at
small angular scale as seen in Fig. 5. Nevertheless, the qualitative features of the cross-correlation functions are unchanged.
Fig. 9 shows the same figure as Fig. 7, but the number densityof the sources is ∼ 10−4 Mpc−3 . The significance and angularscale of positive correlation is similar to the cases of ns ∼ 10−4
Mpc−3 .
So far, we have fixed the maximum distance of HEP sources
considered for correlation as 75 Mpc, which is the same distance as analyses by the PAO (Abraham et al. 2007, 2008) (z = 0.018 in the concordance cosmology ). However, this value
was derived to optimize the correlation signal, and therefore it
does not have inevitability, although it is not inconsistent with
the GZK scenario. Thus, we check the correlation by changing
the maximum distance considered, adopting 100 Mpc and 50
Mpc.
Fig. 10 shows the cross-correlation functions of 200 protons
with energies above 6 × 1019 eV and their sources within 100
Mpc. ns ∼ 10−4 Mpc−3 (upper panels) and ns ∼ 10−5 Mpc−3(lower panels) are assumed. The anisotropic apertures of thePAO (left panels) and TA (right panels) are taken into account.

Compared to Figs. 7 and 9, qualitative features are the same,that is, positive correlation is predicted at small angular scale
except for the AS GMF models in the northern sky. The exception shows the cross-correlation functions almost consistent
with no correlation. Fig. 11 shows the same as Fig. 10, but the
maximum distance of HEP sources considered for correlation
is 50 Mpc. They have also the same qualitative features as Fig.
10. However, we can find that the significance of positive correlation is lower than in Figs. 7, 9 and 10, because of relatively
large error bars. A main reason of this is a smaller volume including considered sources. Since the number of sources used
for correlation analysis is small, the fluctuations of the crosscorrelation functions become large. Thus, the investigation of
spatial correlation with nearby sources is not always efficient.
On the other hand, considering distant sources is also not always better because there are the deflections of HEP trajectories by IGMF, which has much uncertainty, in a real situation.
Therefore, the maximum distance of sources with which thecorrelation is considered should be treated as a parameter atpresent, when we analyze observational data actually.
4. DISCUSSION & CONCLUSION

In this study, we considered the propagation of HEPs in the
Galactic space and investigated the cross-correlation between
the arrival directions of HEPs with energies above 6 × 1019 eV
and the projected positions of their sources by using simulations. We found that the expected correlation between HEPs
and their sources is not disturbed and is positive at small angular scale in many cases after 200 HEPs detection, althougha coherent (spiral) component of GMF efficiently deflects theDoes GMF Disturb the Correlation of HECRs with their Sources?
90o

90o

60o

11
60o

30o

180o

30o

-180o

180o

-180o

-30o

-90o
90o

-30o

-60o

-90o
90o

60o

-60o
60o

30o

180o

30o

-180o

180o

-180o

-30o

-90o

-60o

-30o

-90o

-60o

F IG . 12.— Arrival directions of the PAO events published in Abraham et al. (2008) (red) and the possible directions of their sources on the assumption that all
the PAO events are protons (green) and irons(blue). In other words, we backtrack the PAO events following the BS-S model (upper left), BS-A model (upper right),
AS-S model (lower left), and AS-A model (lower right), and plot their velocity directions at the boundary between the Galactic space and extragalactic space.

trajectories of HEPs. An exception is the cases of the AS GMF
models in the northern sky. In these cases, the predicted crosscorrelation is consistent with no correlation or has very low significance even if it is positive, which reflects that the BS models predict the smaller deflections of HEPs than the AS models
especially in the northern sky. Since the difference of the deflection angles is produced by field reversals of GMF, null observation of the correlation between HEPs and their sources orvery weak correlation imply no field reversals of GMF at leastoutside the solar system. These features are qualitatively unchanged if we change the maximum distance of HEP sources
with which we consider the correlation in the range within 100
Mpc.
Finally we discuss 3 points which are able to affect the results of this study; the compositions of HECRs (possibility(ii) listed in §1), GMF in the Galactic halo, and intergalacticmagnetic field (possibility (iii) listed in §1). Recent reports of
the PAO showed that HECRs include a large fraction of theheavier nuclei by < Xmax > and the deviations of Xmax measurements (Bellido et al. 2009) (although the HiRes reported
proton-dominated composition at the highest energy (Belz etal. 2009)). A remarkable point is that the PAO simultaneously
reported that the correlation between the arrival directions of
HECRs and the positions of nearby astrophysical objects (Abraham et al. 2007, 2008; Hague et al. 2009). Since the trajectories of heavy nuclei are more deflected by intervening magnetic
fields than protons, the observed spatial correlation tends to be
disturbed if heavy nuclei are included in observed HECRs.
It is expected from a reaction with CMB photons that themain component of heavy nuclei in HECRs is iron. Nuclei
can be broken through giant dipole resonance by interactions
with CMB photons at ultra-high energy. The giant dipole resonance predicts a GZK-like spectral steepening, but the energy of steepening is lower than the GZK steepening because
its threshold energy is smaller than that of photopion production of protons (e.g., Allard et al. (2005)). The energy of thesteepening by the giant dipole resonance increases as increas-

ing the atomic number of nuclei, and then the energy becomescomparable with the energy of the GZK steepening. Thus, thedominant component of heavy nuclei in HECRs should be irons
to fit an observed energy spectrum with the GZK-like feature.
Lighter nuclei disturb such a spectral shape. Compositions
other than protons and irons must be a very small fraction. It isquite serious for the spatial correlation observed by the PAO.
In Fig. 12, we demonstrate that source positions of the PAO
events with the IRAS galaxies within z = 0.018 (black dots) assuming that all of the PAO events are protons (green) and irons(blue). We backtrack the detected PAO events (red) in the 4
GMF models and plot their velocity directions just outside the
Galaxy (at 40 kpc from the Galactic center). These circles represent the intrinsic directions of their sources on the assumption of no IGMF. We can obviously find that the inferred directions of iron sources are not consistent with the distribution of
nearby galaxies. It is a fine-tuning problem that the deflection
by IGMF improves the correlation. Thus, pure iron scenario isproblematic in the viewpoint of the spatial correlation. In orderto satisfy the consistency of the results of composition measurements and the spatial correlation observed by the PAO, HECRs
also include significant fraction of protons. The proton component can unveil the positions of HEP sources as we have shownin this study. In this viewpoint, the fraction of irons to protons
is an essential quantity to explain the current result of the PAO.
We consider only spiral fields of GMF in this study. However, vertical components (perpendicular to the Galactic Plane)of GMF have been observed both in the vicinity of the solar system (Han & Qiao 1994) and around the Galactic center (Han2007). A magnetic field in the Galactic halo might also affectthe deflections of HECR trajectories. The existence of Galactic
wind has been indicated by Galactic diffuse soft X-ray observations (Everett et al. 2008). On the other hand, high-sensitivity
observations of several edge-on galaxies with galactic winds
have revealed so-called X-shaped magnetic fields in their halos(Krause 2007; Beck 2009). These facts imply other structured
components of GMF in the halo of our galaxy. This effect to12Takami & Sato

HECR propagation in the Galaxy should be considered and it isa next target of our studies.
In a real situation, HECRs are deflected by not only GMF but
also IGMF. However, our current understanding about IGMF ispoor and IGMF modellings have large uncertainty. A numericalsimulation on cosmological structure formation in local Universe showed a magnetic structure highly concentrated to dense
region (Dolag et al. 2005). In this case, HEPs are almost not
deflected except the directions of clusters of galaxies, since thecross-sectional area occupied by strong IGMF is small. A simple modelling of IGMF by Takami et al. (2006) also predicted
relatively small deflection angles of HEPs and showed that HEP
sources are unveiled by HEPs themselves within a few degree(Takami & Sato 2008a). These 2 IGMF models are based ona realistic matter structure of local Universe constructed by the
IRAS catalog of galaxies. On the other hand, hydrodynamical
simulations on cosmological structure formation predicted relatively strong IGMF in filamentary structures up to 10nG (Siglet al. 2004; Das et al. 2008). A large cross-sectional areaof filamentary structures with ∼ 10 nG fields makes HEPs be
deflected significantly and a large part of information on theirsources is lost. Note that these IGMFs have uncertainty on observer’s positions because their resultant matter distributions
do not correspond to local Universe actually observed. However, an IGMF model of Das et al. (2008) does not lose thecorrelation between the arrival directions of HEPs and largescale structure of nearby Universe (Ryu et al. 2009). More-over, another modelling of IGMF pointed out the possibilityof fake correlation, which means that HECRs are scattered atstrongly magnetized region not including their sources (Kotera& Lemoine 2008). In this case, the correlation with largescale structure may be also conserved. Thus, it is uncertainat present because of the uncertainty of IGMF itself whether
IGMF disturbs possible correlation between the arrival directions of HEPs and the positions of their sources, though thecorrelation with the large-scale structure is expected by many
works. One possibility to know IGMF structure is a projectof Faraday rotation measurement by Square Kilometer Array
3. The Faraday rotation on the whole sky will give us useful
information on magnetic fields in our local Universe.
In this study, we discussed the effects of GMF in cosmic ray
astronomy. As we briefly discussed in this section, the understanding of the heavy components in HECRs and IGMF will bea key to investigate the potential of cosmic ray astronomy.
H.T. thanks to Veniamin Berezinsky for useful comments at
the early stage of this study and Peter Biermann for comments
on GMF in the Galactic halo. This work is supported by World
Premier International Research Center Initiative (WPI Initiative), the Ministry of Education, Culture, Sports, Science andTechnology (MEXT) of Japan. This work is also supported by
Grants-in-Aid for Scientific Research from the MEXT of Japanthrough No.21840019 (H.T.) and No.19104006 (K.S.).

REFERENCES
Abraham, J., et al. [The Pierre Auger Collaboration] 2007, Science, 318, 938
Abraham, J., et al. [The Pierre Auger Collaboration] 2008, Astropart. Phys., 29,
188
Allard, D., Parizot, E., Olinto, A.V., Khan, E., & Goriely, S., 2005, A&A, 443,
L29
Alvarez-Muniz, J., Engel, R., Stanev, T. 2002, ApJ, 572, 185
Beck, R. 2009, Ap&SS, 320, 77
Bellido, J.A. et al. [The Pierre Auger Collaboration], Proc. 31th Int. Cosmic
Ray Conf. (Lodz) #124
Belz, J., et al. [The HiRes Collaboration], Proc. 31th Int. Cosmic Ray Conf.
(Lodz) #854
Bhattacharjee, P., & Sigl, G. 2000, Phys. Rep., 327, 109
Blake, C., Pope, A., Scott, D., & Mobasher, B. 2006, MNRAS, 368, 732
Blasi, P. & de Marco, D. 2004, Astropart. J., 20, 559
Cuoco, A., Hannestad, S., Haugbolle, T., Kachelriess, M., & Serpico, P. D.
2009, ApJ, 702, 825
Das, S., Kang, H., Ryu, D., & Cho, J. 2008, ApJ, 682, 29
Dolag, K., Grasso, D., Springel, V., & Tkachev, I. 2005, J. Cosmology
Astropart. Phys., 01, 009
Everett J. E., Zweibel, W. G., Benjamin, R. A., McCammon, D., Rocks, L., &Gallagher III, J. S. 2008, ApJ, 674, 258
George, M. R., Fabian A. C., Baumgartner, W. H., Mushotzky, R. F., & Tueller,
J. 2008, MNRAS, 388, L59
Ghisellini, G., Ghirlanda, G., Tavecchio, F., Fraternali, F., & Pareschi, G. 2008,
MNRAS, 390, L88
Greisen, K. 1966, Phys. Rev. Lett., 16, 748
Hague, J. et al. [The Pierre Auger Collaboration] 2009, Proc. 31th Int. Cosmic
Ray Conf. (Lodz) # 143
Han, J. L., & Qiao, G.J. 1994, A&A, 288, 759
Han, J.L. 2007, in IAU Symp. 242, Astrophysical Masers and Their
Environments, ed. J.M.Chapman & W.A.Baan (Cambridge: Cambridge
Univ. Press), 55
Harari, D., Mollerach, S., & Roulet, E. 1999, J. High Energy Phys., 08, 022
Kachelriess, M., & Semikoz, D. 2005, Astropart. Phys., 23, 486
Kachelriess, M., Serpico, P. D., Teshima, M. 2007, Astropart. Phys., 26, 378

3 http://www.skatelescope.org/

Kalashev, O. E., Kuz’min, V. A., & Semikoz, D. V. 2001, Mod. Phys. Lett., 16,2505
Kashti, T. & Waxman, E. 2008, J. Cosmology Astropart. Phys., 05, 009
Kotera, K. & Lemoine, M. 2008, Phys. Rev. D, 77, 123003
Krause, M. 2007, Mem. Soc. Astron. Ital., 78, 314
Medina-Tanco, G. A., de Gouveia dal Pino, E. M., & Horvath, J. E. 1998, ApJ,
492, 200
Moskalenko, I. V., Stawarz, L., Porter, T. A., & Cheung, C. C. 2009, ApJ, 693,
1261
Mucke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000,
Com. Phys. Comm., 124, 290
Nonaka, T. et al. [The Telescope Array Collaboration], Proc. 31th Int. Cosmic
Ray Conf. (Lodz) #974
Prouza, M., Smida, R. 2003, A&A, 410, 1
Ryu, D., Kang, H., & Das, S., Proc. 31th Int. Cosmic Ray Conf. (Lodz) #839
Saunders, W. et al. 2000, MNRAS, 317, 55
Sigl, G., Miniati, F., & Ensslin, T.A. 2004, Phys. Rev. D, 70, 043007
Sommers, P. 2001, Astropart. Phys., 14, 271
Stanev, T. 1997, ApJ, 479, 290
Takami, H., Nishimichi, T., Yahata, K., & Sato, K. 2009, J. Cosmology
Astropart. Phys., 06, 031
Takami, H. & Sato, K. 2008a, ApJ, 678, 606
Takami, H. & Sato, K. 2008b, ApJ, 681, 1279
Takami, H., & Sato, K. 2009, Astropart. Phys., 30, 306
Takami, H., Yoshiguchi, H., & Sato, K. 2006, J. Cosmology Astropart. Phys.,639, 803, [Erratum ibid. 653 (2006) 1584]
Takeda, M., et al. [The AGASA Collaboration] 1999, ApJ, 522, 225
Takeuchi, T. T., Yoshikawa, K. & Ishii, T. 2003, ApJ, 587, L89 (erratum ibid.
606 (2004) L171)
Tinyakov, P. G., Tkachev, I. I. 2002, Astropart. Phys., 18, 165
Tinyakov, P. G., & Tkachev, I.I. 2005, Astropart. Phys., 24, 32
Torres, D.& Anchordoqui, L. 2004, Rep. Prog. Phys. 67, 1663
Yoshiguchi, H., Nagataki, S., & Sato, K. 2003, ApJ, 596, 1044
Yoshiguchi, H., Nagataki, S., & Sato, K. 2004, ApJ, 607, 840
Zatsepin, G.T., & Kuz’min, V.A. 1966, JETP Lett., 4, 78