text:
# On astrophysical solution to ultra high energy cosmic rays

## I  Introduction

The systematic study of Ultra High Energy Cosmic Rays (UHECR) started in late 1950s after construction of Volcano Ranch (USA) and Moscow University (USSR) arrays. During the next 50 years of research the origin of UHE particles, which hit the detectors, was not well understood. At present due to the data of the last generation arrays, Haverah Park (UK)M. Ave *et al.* \[Haverah Park collaboration\] (2003), Yakutsk (Russia) V. P. Egorova *et al.* (2004), Akeno and AGASA (Japan) N. Hayashida *et al.* \[AGASA collaboration\] (1994); K. Shinozaki *et al.* \[AGASA collaboration\] (2004), Fly’s Eye D. J. Bird *et al.* \[Fly’s Eye collaboration\] (1994) and HiRes R. U. Abbasi *et al.* \[HiRes collaboration\] (2004) (USA) we are probably very close to understanding the origin of UHECR. The forthcoming data of Auger detector (see The Pierre Auger Collaboration  for the first results) will undoubtedly shed more light on this problem.

On the theoretical side we have an important clue to understanding the UHECR origin: the interaction of extragalactic protons, nuclei and photons with CMB, which leaves the imprint on UHE proton spectrum, most notably in the form of the Greisen-Zatsepin-Kuzmin (GZK) Greisen (1966); Zatsepin and Kuzmin (1966) cutoff for the protons.

We shortly summarize the basic experimental results and the results of the data analysis, important for understanding of UHECR origin (for a review see Nagano and Watson (2000)).
(i) The spectra of UHECR are measured with good accuracy at $`1\times 10^{18}1\times 10^{20}`$ eV, and these data have a power to reject or confirm some models. The discrepancy between the AGASA and HiRes data at $`E>1\times 10^{20}`$ eV might have the statistical origin (see D. De Marco, P. Blasi and Olinto (2003) and discussion in Section IV.3), and the GZK cutoff may exist.
(ii) The mass composition at $`E>1\times 10^{18}`$ eV (as well as below) is not well known (for a review see Watson (2004)).

Different methods result in different mass composition, and the same methods disagree in different experiments. In principle, the most reliable method of measuring the mass composition is given by elongation rate (energy dependence of maximum depth of a shower, $`X_{\mathrm{max}}`$) measured by the fluorescent method. The data of Fly’s Eye in 1994 D. J. Bird *et al.* \[Fly’s Eye collaboration\] (1994) favored iron nuclei at $`1\times 10^{18}`$ eV with a gradual transition to the protons at $`1\times 10^{19}`$ eV. The further development of this method by the HiRes detector, which is the extension of Fly’s Eye, shows the transition to the proton composition already at $`1\times 10^{18}`$ eV Archbold and Sokolsky (2003); R. U. Abbasi *et al.* \[HiRes Collaboration\] (2005). At present the data of HiRes Archbold and Sokolsky (2003); R. U. Abbasi *et al.* \[HiRes Collaboration\] (2005), HiRes-MIA T. Abu-Zaayad *et al.* \[HiRes Collaboration\] (2000) and Yakutsk A. V. Glushkov *et al.* \[Yakutsk collaboration\] (2000) favor the proton-dominant composition at $`E>1\times 10^{18}`$ eV, data of Haverah Park M. Ave *et al.* \[Haverah Park collaboration\] (2003) do not contradict such composition at $`E>(12)\times 10^{18}`$ eV, while data of Fly’s Eye D. J. Bird *et al.* \[Fly’s Eye collaboration\] (1994) and Akeno M. Honda *et al.* \[Akeno Collaboration\] (1993) agree with mixed composition dominated by iron.
(iii) The arrival directions of particles with energy $`E4\times 10^{19}`$ eV show the small-angle clustering within the angular resolution of detectors. AGASA found 3 doublets and one triplet among 47 detected particles M. Takeda \[AGASA collaboration\] (1999) (see the discussion in Finley and Westerhoff (2004)). In the combined data of several arrays Y. Uchihori *et al.* (2000) there were found 8 doublets and 2 triplets in 92 events. The stereo HiRes data R. U. Abbasi *et al.* (2004) do not show small-angle clustering for 27 events at $`E4\times 10^{19}`$ eV, maybe due to limited statistics.

Small-angle clustering is most naturally explained in the case of rectilinear propagation as a random arrival of two (three) particles from a single source S. L. Dubovsky, P. G. Tinyakov and Tkachev (2000). This effect has been calculated in Refs. Fodor and Katz (2001); H. Yoshiguchi, S. Nagataki and Sato (2003); G. Sigl, F. Miniati and Ensslin (2003a); H. Yoshiguchi, S. Nagataki and Sato (2004a); P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a); D. De Marco, P. Blasi and Olinto . In the last five works the calculations have been performed by the Monte Carlo (MC) simulations and results agree well. According to Kachelriess and Semikoz (2005a) the density of the sources, needed to explain the observed number of doublets is $`n_s=(13)\times 10^5`$ Mpc-3. In P. Blasi and D. De Marco (2004) the best fit is given by $`n_s1\times 10^5`$ Mpc-3 and the large uncertainties (in particular due to ones in observational data) are emphasized.
(iv) Recently there have been found the statistically significant correlations between directions of particles with energies $`(48)\times 10^{19}`$ eV and directions to AGN of the special type - BL Lacs Tinyakov and Tkachev (2001) (see also the criticism N. W. Evans, F. Ferrer and Sarkar (2004) and the reply Tinyakov and Tkachev (2004)).

The items (iii) and (iv) favor rectilinear propagation of primaries from the point-like extragalactic sources, presumably AGN. However, the propagation in magnetic fields also exhibits clustering M. Lemoine, G. Sigl and Biermann ; G. Sigl, F. Miniati and Ensslin (2003a); H. Yoshiguchi, S. Nagataki and Sato (2004b).

The quasi-rectilinear propagation of ultra-high energy protons is found possible in MHD simulations K. Dolag, D. Grasso, V. Springel and Tkachev (2004) of magnetic fields in large-scale structures of the universe (see however the simulations in G. Sigl, F. Miniati and Ensslin (2003b); G. Sigl, F. Miniati and Ensslin (2004) with quite different results).

There are many unresolved problems in the field of Ultra High Energy Cosmic Rays, such as nature of primaries (protons? nuclei? or the other particles?), transition from galactic to extragalactic cosmic rays, sources and acceleration, but most intriguing problem remains existence of superGZK particles with energies higher than $`E1\times 10^{20}`$ eV. “The AGASA excess”, namely 11 events with energy higher than $`1\times 10^{20}`$ eV, is still difficult to explain, though there are indications that it may have the statistical origin combined with systematic errors in energy determination (see section IV.3). The AGASA excess, if it is real, should be described by another component of UHECR, most probably connected with the new physics: superheavy dark matter, new signal carriers, like e.g. light stable hadron and strongly interacting neutrino, the Lorentz invariance violation etc.

The problem with superGZK particles is seen in other detectors, too. Apart from the AGASA events, there are five others: the golden Fly’s Eye event with $`E3\times 10^{20}`$ eV, one HiRes event with $`E1.8\times 10^{20}`$ eV and three Yakutsk events with $`E1\times 10^{20}`$ eV. No sources are observed in the direction of these particles at the distance of order of attenuation length. The most severe problem is for the golden Fly’s Eye event: with attenuation length $`l_{\mathrm{att}}=21`$ Mpc and the homogeneous magnetic field 1 nG on this scale, the deflection of particle is only $`3.7^{}`$. Within this angle there are no remarkable sources at distance $`20`$ Mpc Elbert and Sommers (1995).

In this paper we analyze the status of most conservative astrophysical solution to ultra-high energy cosmic ray problem, assuming that primary particles are protons or nuclei accelerated in extragalactic sources. In the first part of the paper (Sections II - V) we analyze the signatures of ultra-high energy protons propagating through CMB. We found that the dip, a spectral feature in energy range $`1\times 10^{18}4\times 10^{19}`$ eV, is well confirmed by observational data of AGASA, HiRes, Yakutsk and Fly’s Eye detectors. In Sections VI -VII we discuss in the model-dependent way the transition from galactic to extragalactic cosmic rays and extragalactic sources: AGN and GRBs.

We use in formulae for the flux throughout the paper, e.g. in equations (10), (11), (13), (14), (20) and in Appendix E, energy $`E`$ measured in GeV, luminosity $`L_p`$ \- in GeV/s, emissivity (comoving energy density production rate) $``$ \- in GeV cm-3 s-1 and $`E_{\mathrm{min}}=1`$ GeV, if not otherwise indicated.

## II  Energy losses and the universal spectrum of UHE protons

In this Section we present our recent calculations of energy losses for UHE protons interacting with CMB, and calculate the spectrum of protons, assuming the homogeneous distribution of sources in the space and continuous energy losses. The spectrum is calculated, using the conservation of the number of protons, interacting with CMB. Formally it does not depend on the mode of proton propagation (e.g. rectilinear or diffusive), and we shall discuss when this approximation is valid. The proton spectrum calculated in this way we call universal.

### II.1 Energy losses

We present here the accurate calculations of energy losses due to pair production, $`p+\gamma _{CMB}p+e^++e^{}`$, and pion production, $`p+\gamma _{CMB}N+\mathrm{pions}`$, where $`\gamma _{CMB}`$ is a microwave photon.

The energy losses of UHE proton per unit time due to its interaction with low energy photons is given by

$$\frac{1}{E}\frac{dE}{dt}=\frac{c}{2\mathrm{\Gamma }^2}_{\epsilon _{th}}^{\mathrm{}}𝑑\epsilon _r\sigma (\epsilon _r)f(\epsilon _r)\epsilon _r_{\epsilon _r/2\mathrm{\Gamma }}^{\mathrm{}}𝑑\epsilon \frac{n_\gamma (\epsilon )}{\epsilon ^2},$$
(1)
where $`\mathrm{\Gamma }`$ is the Lorentz factor of the proton, $`\epsilon _r`$ is the energy of background photon in the reference system of the proton at rest, $`\epsilon _{th}`$ is the threshold of the considered reaction in the rest system of the proton, $`\sigma (\epsilon _r)`$ is the cross-section, $`f(\epsilon _r)`$ is the mean fraction of energy lost by the proton in one $`p\gamma `$ collision in the laboratory system, (see Fig. 22 in Appendix A), $`\epsilon `$ is the energy of the background photon in the laboratory system, and $`n_\gamma (\epsilon )`$ is the density of background photons.

For CMB with temperature $`T`$ Eq. (1) is simplified

$$\frac{1}{E}\frac{dE}{dt}=\frac{cT}{2\pi ^2\mathrm{\Gamma }^2}_{\epsilon _{th}}^{\mathrm{}}𝑑\epsilon _r\sigma (\epsilon _r)f(\epsilon _r)\epsilon _r\left\{\mathrm{ln}\left[1\mathrm{exp}\left(\frac{\epsilon _r}{2\mathrm{\Gamma }T}\right)\right]\right\}.$$
(2)
Further on we shall use the notation

$$\beta _0(E)=\frac{1}{E}\frac{dE}{dt},b_0(E)=\frac{dE}{dt}$$
(3)
for energy losses on CMB at present cosmological epoch, $`z=0`$ and $`T=T_0`$. For the epochs with red-shift $`z`$ one has:

$`\beta (E,z)`$ $`=`$ $`(1+z)^3\beta _0[(1+z)E],`$ (4)
$`b(E,z)`$ $`=`$ $`(1+z)^2b_0[(1+z)E].`$ (5)
Another important physical quantity needed for calculations of spectra is the derivative $`db_0(E)/dE`$, which can be calculated as

$$\frac{db_0(E)}{dE}=\beta _0(E)+\frac{c}{4\pi ^2\mathrm{\Gamma }^3}_{\epsilon _{\mathrm{th}}}^{\mathrm{}}𝑑\epsilon _r\sigma (\epsilon _r)f(\epsilon _r)\frac{\epsilon _r^2}{\mathrm{exp}\left(\frac{\epsilon _r}{2\mathrm{\Gamma }T_0}\right)1}.$$
(6)
As one can see from Fig. 1, $`db_0(E)/dE`$ is numerically very similar to $`\beta _0(E)`$, and for approximate calculations one can use $`\beta _0(E)`$ values for both functions.

From Eqs. (1) and (2) one can see that apart from cross-section the mean fraction of energy lost by the proton in laboratory system in one collision, $`f(\epsilon _r)`$ (see Eq. (28)), is the basic quantity needed for calculations of energy losses. The threshold values of these quantities are well known:

$$f_{\mathrm{pair}}=\frac{2m_e}{2m_e+m_p},f_{\mathrm{pion}}=\frac{m_\pi }{m_\pi +m_p},$$
(7)
where $`f_{\mathrm{pair}}`$ and $`f_{\mathrm{pion}}`$ are the threshold fractions for $`p+\gamma p+e^++e^{}`$ and $`p+\gamma N+\pi `$, respectively.

Pair production loss has been previously discussed in many papers. All authors directly or indirectly followed the standard approach of Ref. Blumenthal (1970) where the first Born approximation of the Bethe-Heitler cross-section with proton mass $`m_p\mathrm{}`$ was used. In contrast to Ref. Blumenthal (1970), we use here the first Born approximation approach of Ref. Berg and Linder (1961) accounting for the finite proton mass. This allows us to calculate the average fraction of energy lost by the proton in laboratory system by performing a fourfold integration over invariant mass of electron-positron pair, $`M_X`$, over an angle between incident and scattered proton, and polar and azimuthal angles of an electron in the c.m.s. of the $`e^+e^{}`$-pair (see Appendix A for further details).

Calculating photopion energy loss we follow methods developed in papers Berezinsky and Gazizov (1993); Gazizov (1996). Total cross-sections are taken according to Ref. Gabathuler (1974). At low c.m.s. energy, $`E_c`$, we consider the binary reactions $`p+\gamma \pi ^0+p`$, $`p+\gamma \pi ^++n`$ (they include the resonance $`p+\gamma \mathrm{\Delta }^+`$). Differential cross-sections of binary processes at small energies are taken from D. Menze, W. Pfeil and Wilke (1977); Wolf and Söding (1978); Fisher (1974). At $`E_c>4.3`$ GeV we assume the scaling behavior of differential cross-sections, the latter being taken from Ref. Meyer (1974); Brasse (1974); K. C. Moffeit *et al.* (1972). In the intermediate energy range we interpolate between the angular distribution of these two regimes with the cross-section being the difference between total cross-section and cross-section of the all binary processes. Angular distributions for this part of cross-section vary from isotropic at threshold to those imposed by inclusive pion photoproduction data at high energies. The overall differential cross-sections coincide with low-energy binary description and high-energy scaling distributions and join smoothly these two regimes in the intermediate region.

The results of our calculations are presented in Fig. 1 in terms of $`E^1dE/dt`$ as function of energy (curve 1). Also plotted is the derivative $`db_0(E)/dE`$ (Fig. 1b). This quantity is needed for calculation of differential energy spectrum. In Fig. 1 we plot for comparison the energy losses as calculated by Berezinsky and Grigorieva 1988 Berezinsky and Grigorieva (1988) (dashed curve 2). The difference in energy losses due to pion production is very small, not exceeding 5% in the energy region relevant for comparison with experimental data ($`E10^{21}`$eV). The difference with energy losses due to pair production is larger and reaches maximal value 15%. The results of calculations by Stanev *et al.* T. Stanev, R. Engel, A. Muecke, R. J. Protheroe and Rachen (2000) are shown by black squares. These authors have performed the detailed calculations for both aforementioned processes, though their approach is different from ours, especially for photopion process. Our new energy losses are practically indistinguishable from Stanev *et al.* T. Stanev, R. Engel, A. Muecke, R. J. Protheroe and Rachen (2000) for pair production and pion production at low energies, and differ by 15-20% for pion production at highest energies (see Fig. 1). Stanev *et al.* claimed that energy losses due to pair production is underestimated by Berezinsky and Grigorieva Berezinsky and Grigorieva (1988) by 30-40%. Comparison of data files of Stanev *et al.* and Berezinsky and Grigorieva (see also Fig. 1) shows that this difference is significantly less. Most probably, Stanev *et al.* scanned inaccurately the data from the journal version of the paper Berezinsky and Grigorieva (1988).

### II.2 Universal spectrum

To calculate the spectrum one should first of all evolve the proton energy $`E`$ from the time of observation $`t=t_0`$ (or $`z=0`$) to the cosmological epoch of generation $`t`$ (or red-shift $`z`$), using the the adiabatic energy losses $`EH(z)`$ and $`b(E,z)`$ given by Eq. (4):

$$dE/dt=EH(z)+(1+z)^2b_0\left[(1+z)E\right],$$
(8)
where $`H(z)=H_0\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}`$ is the Hubble parameter at cosmological epoch $`z`$, with $`H_0=72`$ km/s Mpc, $`\mathrm{\Omega }_m=0.27`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.73`$ D. N. Spergel *et al.* \[WMAP Collaboration\] (2003).

We calculate the spectrum from conservation of number of particles in the comoving volume (protons change their energy but do not disappear). For the number of UHE protons per unit comoving volume, $`n_p(E)`$, one has:

$$n_p(E)dE=_{t_{\mathrm{min}}}^{t_0}𝑑tQ_{\mathrm{gen}}(E_g,t)𝑑E_g,$$
(9)
where $`t`$ is an age of the universe, $`E_g=E_g(E,t)`$ is a generation energy at age $`t`$, calculated according to Eq. (8) and $`Q_{\mathrm{gen}}(E_g,t)`$ is the generation rate per unit comoving volume, which can be expressed through emissivity $`_0`$, the energy release per unit time and unit of comoving volume, at $`t=t_0`$, as

$$Q_{\mathrm{gen}}(E_g,t)=_0(1+z)^mKq_{\mathrm{gen}}(E_g),$$
(10)
where $`(1+z)^m`$ describes the possible cosmological evolution of the sources. In the case of the power-law generation, $`q_{\mathrm{gen}}(E_g)=E_g^{\gamma _g}`$, with normalization constant $`K=\gamma _g2`$ for $`\gamma _g>2`$ and $`K=(\mathrm{ln}E_{\mathrm{max}}/E_{\mathrm{min}})^1`$ for $`\gamma _g=2`$. We recall that in these formulae and everywhere below energies $`E`$ are given in GeV, emissivity $``$ in GeV cm-3s-1 and source luminosity $`L`$ in GeV s-1.

From Eq. (9) one obtains the diffuse flux as

$$J_p(E)=\frac{c}{4\pi }_0K\times _0^{z_{max}}𝑑z\left|\frac{dt}{dz}\right|(1+z)^mq_{\mathrm{gen}}(E_g)\frac{dE_g}{dE},$$
(11)
where $`|dt/dz|=H^1(z)/(1+z)`$ and analytic expression for $`dE_g/dE`$ is given by Eq. (36) in Appendix B.

The spectrum (11) is referred to as universal spectrum. Formally it is derived from conservation of number of particles and does not depend on propagation mode (see Eq. (9)). But in fact, the homogeneity of the particles, tacitly assumed in this derivation, implies the homogeneity of the sources, and thus the condition of validity of universal spectrum is a small separation of sources. The homogeneous distribution of particles in case of homogeneous distribution of sources and inhomogeneous magnetic fields follows from the Liouville theorem (see Ref. Aloisio and Berezinsky (2004)).

Several effects could in principle modify the shape of universal spectrum. They include propagation in magnetic fields, discreteness in distribution of the sources, large-scale inhomogeneous distribution of sources, local source overdensity or deficit, and fluctuations in interaction. These effects will be studied in the next sections. With exception of energies beyond GZK cutoff and energies below $`1\times 10^{18}`$ eV, the universal spectrum is only weakly modified by aforementioned effects. Here we shortly comment on role of magnetic fields.

As numerical simulations (see e.g. G. Sigl, M. Lemoine and Biermann (1999); H. Yoshiguchi *et al.* (2003)) show, the propagation of UHE protons in strong magnetic fields changes the energy spectrum (for physical explanation of this effect see Aloisio and Berezinsky (2004)). The influence of magnetic field on spectrum depends on the separation of the sources, $`d`$. For uniform distribution of sources with separation $`d`$ much less than characteristic lengths of propagation, such as attenuation length $`l_{\mathrm{att}}`$ and the diffusion length $`l_{\mathrm{diff}}`$, the diffuse spectrum of UHECR is the universal one independent of mode of propagation Aloisio and Berezinsky (2004). This statement has a status of the theorem. For the wide range of magnetic fields 0.1 - 10 nG and distances between sources $`d<50`$ Mpc the spectrum at $`E>10^{18}`$ eV is close to the universal one Aloisio and Berezinsky (2005).

### II.3 Modification factor

The analysis of spectra is very convenient to perform in terms of the modification factor.
Modification factor is defined as a ratio of the spectrum $`J_p(E)`$ (see Eq. (11)) with all energy losses taken into account, to the unmodified spectrum $`J_p^{\mathrm{unm}}(E)`$, where only adiabatic energy losses (red shift) are included.

$$\eta (E)=\frac{J_p(E)}{J_p^{\mathrm{unm}}(E)}.$$
(12)
For the power-law generation spectrum for $`\gamma _g>2`$ without evolution one has

$$J_p^{\mathrm{unm}}(E)=\frac{c}{4\pi }(\gamma _g2)_0E^{\gamma _g}_0^{z_{\mathrm{max}}}𝑑z\left|\frac{dt}{dz}\right|(1+z)^{\gamma _g+1}.$$
(13)
The modification factor is a less model-dependent quantity than the spectrum. In particular, it should depend weakly on $`\gamma _g`$, because both numerator and denominator in Eq. (12) include $`E^{\gamma _g}`$. In the next Section we consider the non-evolutionary case $`m=0`$ (see Section IV.5 for discussion of evolution). In Fig. 2 the modification factor is shown as a function of energy for two spectrum indices $`\gamma _g=2.0`$ and $`\gamma _g=2.7`$. As expected above, they do not differ much from each other. Note that by definition $`\eta (E)1`$.

## III  Signatures of UHE protons interacting with CMB

The extragalactic protons propagating through CMB have signatures in the form of three spectrum features: GZK cutoff, dip and bump. The dip is produced due to $`e^+e^{}`$-production and bump – by pile-up protons accumulated near beginning of the GZK cutoff. We add here the fourth signature: the second dip.

The analysis of these features, especially dip and bump, is convenient to perform in terms of modification factor Berezinsky and Grigorieva (1988); T. Stanev *et al.* (2000). For the GZK cutoff we shall use the traditional spectra.

### III.1 GZK cutoff

The GZK cutoff Greisen (1966); Zatsepin and Kuzmin (1966) is most remarkable phenomenon, which describes the sharp steepening of the spectrum due to pion production. The GZK cutoff is a model–dependent feature of the spectrum, e.g. the GZK cutoff for a single source depends on the distance to the source.

A common convention is that the GZK cutoff is defined for diffuse flux from the sources uniformly distributed over the universe. In this case one can give two definitions of the GZK-cutoff position. In the first one it is determined as the energy, $`E_{GZK}4\times 10^{19}eV`$, where the steep increase in the energy losses starts (see Fig. 1). The GZK cutoff starts at this energy. The corresponding path length of a proton is $`R_{GZK}(E^1dE/cdt)^11.3\times 10^3`$ Mpc. The advantage of this definition of the cutoff energy is independence of a spectrum index, but this energy is too low to judge about presence or absence of the cutoff in the measured spectrum. More practical definition is $`E_{1/2}`$, where the flux with cutoff becomes lower by factor 2 than power-law extrapolation. This definition is convenient to use for the integral spectrum, which is better approximated by power-law function, than the differential one.

In Fig. 3 the function $`E^{(\gamma 1)}J(>E)`$, where $`J(>E)`$, the calculated integral diffuse spectrum, is plotted as function of energy. Note, that $`\gamma >\gamma _g`$ is an effective index of the power-law approximation of the spectrum modified by energy losses. For wide range of generation indices $`2.1\gamma _g2.7`$ the cutoff energy is the same, $`E_{1/2}5.3\times 10^{19}`$ eV. The corresponding proton path length is $`R_{1/2}800`$ Mpc. In panel a) of Fig. 4 $`E_{1/2}`$ is found from the integral spectrum of the Yakutsk array in the reasonable agreement with theoretical prediction. The HiRes data are shown in panel b). These data have large uncertainties, which prevent the accurate determination of $`E_{1/2}`$. However, they agree with the predicted value $`E_{1/2}5.3\times 10^{19}`$ eV. In the recent paper D. R. Bergman *et al.* \[HiRes Collaboration\]  the Hires collaboration found $`E_{1/2}=(5.9_{0.8}^{+2.4})\times 10^{19}`$ eV in better agreement with the predicted value.

In Fig. 5 the calculated universal spectra with the GZK cutoff are compared with AGASA and HiRes data N. Hayashida *et al.* \[AGASA collaboration\] (1994); K. Shinozaki *et al.* \[AGASA collaboration\] (2004); R. U. Abbasi *et al.* \[HiRes collaboration\] (2004). While the HiRes data agree with the predicted GZK cutoff, the AGASA data show significant excess over this prediction at $`E>1\times 10^{20}`$ eV.

The GZK cutoff as calculated above has many uncertainties.

The energy shape of the GZK feature is model dependent. The local excess of sources makes it flatter, and the deficit – steeper. The shape is affected by $`E_{\mathrm{max}}^{\mathrm{acc}}`$ (see Fig. 5) and by fluctuations of source luminosities and distances between the sources. The cutoff, if discovered, can be produced as the acceleration cutoff (steepening below the maximum energy of acceleration in the generation spectrum). Since the shape of both, the GZK cutoff and acceleration cutoff, is model-dependent, it will be difficult to argue in favor of any of them, in case a cutoff is discovered.

To illustrate the effect of local overdensity/deficit of the UHECR sources we calculate here the UHECR spectra with different local ratios $`n/n_0`$, where $`n`$ is the local density of the UHECR sources and $`n_0`$ is mean extragalactic source density. We use the various sizes of overdensity/deficit regions $`R_{overd}`$, equal to 10, 20 and 30 Mpc. The results of our calculations for the overdensity case are presented in Fig. 6 for $`\gamma _g=2.7`$, $`m=0`$ and for four values of overdensity $`n/n_0`$ equal to 1, 2, 3 and 10, assuming the size of overdensity region 30 Mpc (the results for $`R_{overd}=20`$Mpc are not much different). The spectra for the case of the deficit are shown by curves 5 and 6 in the both panels. The theoretical spectra shown in Fig. 6 illustrate uncertainties in the prediction of the shape of the GZK feature.

An interesting question is whether the local overdensity of the sources can explain the AGASA excess. This problem has been already addressed in M. Blanton, P. Blasi and Olinto (2001) for the realistic distribution of visible galaxies, considering them as the UHECR sources. Our calculations show (see left panel of Fig. 6) that unrealistic overdensity $`n/n_010`$ is needed to explain the AGASA excess.

Both effects, $`E_{\mathrm{max}}`$ and local source overdensity (deficit) affect weakly the shape of the GZK cutoff at $`E<1\times 10^{20}`$ eV. Thus, the precise measurements of the spectrum in this energy region, as well as measurement of $`E_{1/2}`$, give the best test of the GZK cutoff. At higher energies theoretical predictions have large model-dependent uncertainties.

### III.2 Bump in the diffuse spectrum

Protons do not disappear in the photopion interactions, they only loose energy and are accumulated near beginning of the GZK cutoff in the form of a bump.

We see no indication of the bump in the modification-factor energy dependence in Fig. 2. As explained above, it should have been located at merging of $`\eta _{ee}(E)`$ and $`\eta _{\mathrm{tot}}(E)`$ curves. The absence of the bump in the diffuse spectrum can be easily understood. The bumps are clearly seen in the spectra of the single remote sources (see left panel in Fig. 7). These bumps, located at different energies, produce a flat feature, when they are summed up in the diffuse spectrum. This effect is illustrated by Fig. 7 (right panel). The diffuse flux there is calculated in the model where sources are distributed uniformly in the sphere of radius $`R_{\mathrm{max}}`$ (or $`z_{\mathrm{max}}`$). When $`z_{\mathrm{max}}`$ are small (between 0.01 and 0.1) the bumps are seen in the diffuse spectra. When radius of the sphere becomes larger, the bumps merge producing the flat feature in the spectrum. If the diffuse spectrum is plotted as $`E^3J_p(E)`$ this flat feature looks like a pseudo-bump.

### III.3 The dip

The dip is more reliable signature of interaction of protons with CMB than GZK feature. The shape of the GZK feature is strongly model-dependent (see Section III.1), while the shape of the dip is fixed and has a specific form which is difficult to mimic by other mechanisms, unless they have many free parameters. The protons in the dip are collected from the large volume with the radius about 1000 Mpc and therefore the assumption of the uniform distribution of sources within this volume is well justified, in contrast to the GZK cutoff, which strongly depends on local overdensity or deficit of the sources. The GZK cutoff can be mimicked by acceleration cutoff, and since the shape of the GZK cutoff is not reliably predicted, these two cases could be difficult to distinguish.

The problem of identification of the dip depends on the accuracy of observational data, which should confirm the specific (and well predicted) shape of this feature. Do the present data have the needed accuracy?

The comparison of the calculated modification factor with that obtained from the Akeno-AGASA data, using $`\gamma _g=2.7`$, is given in Fig. 8. It shows the excellent agreement between predicted and observed modification factors for the dip.

In Fig. 8 one observes that at $`E<1\times 10^{18}`$ eV the agreement between calculated and observed modification factors becomes worse and the observational modification factor becomes larger than 1. Since by definition $`\eta (E)1`$, it signals the appearance of another component of cosmic rays, which is almost undoubtedly given by galactic cosmic rays. The condition $`\eta >1`$ implies the dominance of the new (galactic) component, the transition occurs at $`E<1\times 10^{18}`$ eV.

To calculate $`\chi ^2`$ for the confirmation of the dip by Akeno-AGASA data, we choose the energy interval between $`1\times 10^{18}`$ eV and $`4\times 10^{19}`$ eV (the energy of intersection of $`\eta _{ee}(E)`$ and $`\eta _{\mathrm{tot}}(E)`$). In calculations we used the Gaussian statistics for low-energy bins, and the Poisson statistics for the high energy bins of AGASA. It results in $`\chi ^2=19.06`$. The number of Akeno-AGASA bins is 19. We use in calculations two free parameters: $`\gamma _g`$ and the total normalization of spectrum. In effect, the confirmation of the dip is characterized by $`\chi ^2=19.06`$ for d.o.f.=17, or $`\chi ^2`$/d.o.f.=1.12, very close to the ideal value 1.0 for the Poisson statistics.

In the right-upper panel of Fig. 8 the comparison of modification factor with the HiRes data is shown. The agreement is also very good: $`\chi ^2=19.5`$ for $`d.o.f.=19`$ for the Poisson statistics. The Yakutsk and Fly’s eye data (not shown here) agree with dip as well. The Auger spectrum The Pierre Auger Collaboration  at this preliminary stage does not contradict the dip.

The good agreement of the shape of the dip $`\eta _{ee}(E)`$ with observations is a strong evidence for extragalactic protons interacting with CMB. This evidence is confirmed by the HiRes data on the mass composition Archbold and Sokolsky (2003); R. U. Abbasi *et al.* \[HiRes Collaboration\] (2005). While the data of the Yakutsk array A. V. Glushkov *et al.* \[Yakutsk collaboration\] (2000) and HiRes-MIA T. Abu-Zaayad *et al.* \[HiRes Collaboration\] (2000) support this mass composition, and Haverah Park data M. Ave *et al.* \[Haverah Park collaboration\] (2003) do not contradict it at $`E(12)\times 10^{18}`$ eV, the data of Akeno M. Honda *et al.* \[Akeno Collaboration\] (1993) and Fly’s Eye D. J. Bird *et al.* \[Fly’s Eye collaboration\] (1994) favor the mixed composition dominated by heavy nuclei.

The observation of the dip should be considered as independent evidence in favor of proton-dominated primary composition in the energy range $`1\times 10^{18}4\times 10^{19}`$ eV.

### III.4 The second dip

The second dip in the spectrum of extragalactic UHE protons appears at energy $`E=6.3\times 10^{19}`$ due to interplay between pair production and photopion production. It is the direct consequence of energy $`E_{\mathrm{eq2}}=6.05\times 10^{19}`$ eV, where energy losses due to $`e^+e^{}`$-production and pion-production become equal (see Fig. 1). This spectrum feature is explained as follows.

The pion-production energy loss increases with energy very fast, and at energy slightly below $`E_{\mathrm{eq2}}`$ $`e^+e^{}`$-production dominates and spectrum can be with high accuracy described in continuous energy-loss approximation. At energy slightly higher than $`E_{\mathrm{eq2}}`$ the pion-production dominates and the precise calculation of spectrum should be performed in the kinetic-equation approach. In this method the evolution of number of particles in interval $`dE`$ is given by two compensating terms, describing the particle exit and regeneration due to $`p\gamma `$ collisions. The small continuous energy losses affect only the exit term and break this compensation, diminishing the flux. The exact calculations are given in Appendix D. The second dip is very narrow and its amplitude at maximum reaches $`10\%`$ (see Fig. 24). This feature can be observed by detectors with very good energy resolution, and it gives the precise mark for energy calibration of a detector. It can be observed only marginally by the Auger detector.

## IV  Robustness of the dip prediction

We calculated the dip for the universal spectrum, i.e. for the case when distances between sources are small enough, and the spectrum does not depend on the propagation mode. In this Section we shall study stability of the dip relative to other possible phenomena, namely, discreteness in the source distribution, propagation in magnetic fields etc. We shall consider also some phenomena related to existence of the dip.

### IV.1 Discreteness in the source distribution

As it follows from analysis of the small-scale anisotropy (see S. L. Dubovsky, P. G. Tinyakov and Tkachev (2000); Fodor and Katz (2001); H. Yoshiguchi, S. Nagataki and Sato (2003); G. Sigl, F. Miniati and Ensslin (2003a); H. Yoshiguchi, S. Nagataki and Sato (2004a)), the average distance between UHECR sources is $`d3050`$ Mpc P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a). Such discreteness affects the spectrum, especially at highest energies, when energy attenuation length is comparable with $`d`$.

In this subsection we demonstrate the stability of the dip relative to discreteness of the sources. We illustrate the effect of discreteness by an example of UHE protons propagating rectilinearly from sources located in the vertices of a 3D cubic lattice with spacing $`d`$. Positions of sources are given by coordinates $`x=id`$, $`y=jd`$ and $`z=kd`$, where $`i,j,k=0,\pm 1,\pm 2\mathrm{}`$. The observer is assumed at $`x=y=z=0`$ with no source there. The diffuse flux for the power-law generation spectrum $`E_g^{\gamma _g}`$ is given by summation over all vertices. The maximum distance is defined by the maximum red-shift $`z_{max}`$. Then the observed flux is given by

$$J_p(E)=\frac{(\gamma _g2)_0d}{(4\pi )^2}\underset{i,j,k}{}\frac{E_g(E,z_{ijk})^{\gamma _g}}{(i^2+j^2+k^2)(1+z_{ijk})}\frac{dE_g}{dE},$$
(14)
where $`_0=L_p/d^3`$ is emissivity, $`z_{ijk}`$ is the red-shift for a source with coordinates $`i,j,k`$, and factor $`(1+z_{ijk})`$ takes into account the time dilation.

The calculated spectra for $`d=1,5,10,20,40`$ and $`60`$ Mpc are shown in Fig. 9 in comparison with the AGASA-Akeno data. In calculations we used $`E_{\mathrm{max}}=1\times 10^{22}`$ eV, $`m=0`$ (no evolution), $`z_{\mathrm{max}}=4`$ and $`\gamma _g=2.7`$. Emissivity $`_0`$ is chosen to fit the AGASA data. One can see that discreteness in the source distribution affects weakly the dip, but the effect is more noticeable for the shape of the GZK cutoff.

With $`d`$ decreasing, the calculated spectra regularly converge to the universal one, as it should be according to propagation theorem Aloisio and Berezinsky (2004). This theorem ensures also that the spectra from Fig. 14 are valid for the case of weak magnetic field, when the diffusion length $`l_{\mathrm{diff}}>d`$.

### IV.2 Dip in the case of diffusive propagation

The dip, seen in the universal spectrum, is also present in the case of diffusive propagation in magnetic field Aloisio and Berezinsky (2005). The calculations are performed for diffusion in random magnetic fields with the coherent magnetic field $`B_0=1`$ nG

and up to $`B_0=100`$ nG on the basic scale $`l_c=1`$ Mpc. The calculated spectrum is shown in Fig. 10 for the case $`(B_0,l_c)=(1\mathrm{nG},1\mathrm{Mpc})`$ and distance between sources $`d=50`$ Mpc. The critical energy, where diffusion changes its regime is $`E_b1\times 10^{18}`$ eV. The spectra are shown for three different regimes at $`E<E_b`$ as indicated in Fig. 10. The universal spectrum is also presented. One can see that the dip agrees well with universal spectrum and observational data, while the shape of the GZK cutoff differs considerably from the universal spectrum.

These calculations demonstrate stability of the dip relative to changing of the propagation mode, and sensitivity of the GZK cutoff to the way of propagation.

### IV.3 Energy calibration of the detectors using the dip

Since the position and shape of the dip is robustly fixed by proton interaction with CMB, it can be used as energy calibrator for the detectors. We use the following procedure for the calibration. Assuming the energy-independent systematic error, we shift the energies in each given experiment by factor $`\lambda `$ to reach minimum $`\chi ^2`$ for comparison with the calculated dip. The systematic errors in energy determination of existing detectors exceed 20%, and it determines the expected value of $`\lambda `$. The described procedure results in $`\lambda _A=0.9`$, $`\lambda _{\mathrm{Ya}}=0.75`$ and $`\lambda _{\mathrm{Hi}}=1.2`$ for AGASA, Yakutsk and HiRes detectors, respectively. After this energy calibration the fluxes in all experiments agree with each other.

First we consider the AGASA and HiRes data. There are two discrepancies between these data (see the upper-left panel of Fig. 11): one is described by factor 1.5 - 2.0 in energy region $`1\times 10^{18}8\times 10^{19}`$ eV, and the second – at $`E1\times 10^{20}`$ eV. In Fig. 11 the spectra of Akeno-AGASA and HiRes are shown before and after the energy calibration. One can see the good agreement of the calibrated data at $`E<1\times 10^{20}`$ eV and their consistency at $`E>1\times 10^{20}`$ eV. This result should be considered together with calculations D. De Marco, P. Blasi and Olinto (2006), where it was demonstrated that 11 superGZK AGASA events can be simulated by the spectrum with GZK cutoff in case of 30% error in energy determination. We may tentatively conclude that existing discrepancy between AGASA and HiRes spectra at all energies are due to systematic energy errors and statistics.

In Fig. 11 the energy spectra are shown for AGASA and Yakutsk spectra before and after energy calibration. Again, the best fit to the dip shape results in excellent agreement in the absolute values of fluxes.

The agreement between spectra of all three detectors after energy calibration by the dip confirms the dip as the spectrum feature produced by interaction of the protons with CMB, and demonstrates compatibility of fluxes measured by AGASA, HiRes and Yakutsk detectors.

### IV.4 Dip and extragalactic UHE nuclei

The proton dip has very good agreement with observations (see Fig. 8). However, in all astrophysical sources the nuclei must be also accelerated to the energies, naively, Z times higher than that for protons. Do UHE nuclei in primary flux upset the good agreement seen in Fig. 8? This problem has been recently considered in V. Berezinsky, A. Gazizov and Grigorieva (2005); D. Allard, E. Parizot and Olinto (a); Sigl and Armengaud (2005) (for study of propagation UHE nuclei through CMB see e.g. F. W. Stecker and Salamon (1992); Sigl (2004); T. Yamamoto *et al.* (2004)).

The presence of nuclei in primary extragalactic spectrum modifies the dip V. Berezinsky, A. Gazizov and Grigorieva (2005). In Fig. 12 the dips for helium and iron nuclei are shown in comparison with that for protons. From this figure one can see that presence of 15 - 20 % of nuclei in the primary flux breaks the good agreement of proton deep with observations. The modification factor for cosmic rays composed of protons and nuclei with the fraction $`\lambda _A=Q_A(E)/Q_p(E)`$, where $`Q(E)`$ is generation function for nuclei (A) and protons (p), can be easily calculated as

$$\eta _{\mathrm{tot}}(E)=\frac{\eta _p(E)+\lambda \eta _A(E)}{1+\lambda }.$$
(15)
The fraction $`\lambda `$ is a model-dependent value, which depends on ratio of number densities of gas components $`n_A/n_H`$ in media, where acceleration operates, on ionization of the gas and on the injection mechanism of acceleration V. Berezinsky . Besides, the UHE nuclei can be destroyed inside the source or in its vicinity Sigl and Armengaud (2005).

The strongest distortion of proton modification factor is given by helium nuclei, for which $`n_{\mathrm{He}}/n_\mathrm{H}0.08`$ corresponding to the helium mass fraction $`Y_p=0.24`$. In Fig. 13 the modification factors are shown for the mixed composition of protons and helium with the mixing parameter $`\lambda =0.1`$ (left panel) and $`\lambda =0.2`$ (right panel). One can judge from these graphs about allowed values of mixing parameter $`\lambda `$. If agreement of the proton dip with observations is not incidental (the probability of this is small according to small $`\chi ^2`$/d.o.f.), Fig. 13 should be interpreted as indication to possible acceleration mechanism V. Berezinsky .

### IV.5 Dip and cosmological evolution of the sources

The cosmological evolution of the sources, i.e. increase of the luminosities and/or space densities with red-shift $`z`$, is observed for many astronomical populations. The evolution is reliably observed for star formation rate in the normal galaxies, but this case is irrelevant for UHECR, because neither stars nor normal galaxies can be the UHECR sources due to low cosmic-ray luminosities $`L_p`$ and maximum energy of acceleration $`E_{\mathrm{max}}`$.

AGN, which satisfy these requirements, also exhibit the evolution seen in radio, optical and X-ray observations. The X-ray radiation is probably most relevant tracer for evolution of UHECR because both radiations are feed by the energy release provided by accretion to a massive black hole: X-rays – through radiation of accretion disk, and UHECR – through acceleration in the jets. According to recent detailed analysis in Ueda et al. (2003) and A. J. Barger *et al.* (2005) the evolution of AGN seen in X ray radiation can be described by factor $`(1+z)^m`$ up to $`z_c1.2`$ and is saturated at larger $`z`$. In Ueda et al. (2003) the pure luminosity evolution and pure density evolution are allowed with $`m=2.7`$ and $`m=4.2`$, respectively, and with $`z_c1.2`$ for both cases. In A. J. Barger *et al.* (2005) the pure luminosity evolution is considered as preferable with $`m=3.2`$ and $`z_c=1.2`$. These authors do not distinguish between different morphological types of AGN. It is possible that some AGN undergo weak cosmological evolution, or no evolution at all. BL Lacs Tinyakov and Tkachev (2001), which are important as potential UHECR sources, show no signs of positive cosmological evolution S. L. Morris *et al.* (2003).

In case of UHECR there is no need to distinguish between luminosity and density evolution, because the diffuse flux is determined by the emissivity, which includes both luminosity and space evolution, as it follows from Eq. (10).

In Fig. 14 we present the calculated dip spectrum in evolutionary models, inspired by the data cited above. For comparison we show also the case of absence of evolution $`m=0`$, as can be valid for BL Lacs. From Fig. 14 one can see that the spectra with evolution up to $`z_c>1`$ can explain the observational data down to (a few)$`\times 10^{17}`$ eV and even below, in accordance with early calculations V. Berezinsky, A. Z. Gazizov and Grigorieva (a); Scully and Stecker (2002); Z. Fodor, S. D. Katz, A. Ringwald and Tu (2003) (see D. R. Bergman *et al.*  for recent analysis). However, for any reasonable magnetic fields protons with these energies have small diffusion lengths and the spectrum acquires the diffusion ’cutoff’ at energy $`E_b=1\times 10^{18}`$ eV (see Section IV.2).

We conclude that for many reasonable evolution regimes the dip agrees with observational data as well as the non-evolutionary case $`m=0`$.

## V  Role of interaction fluctuations

UHE proton spectrum is affected by fluctuations in the photopion production. These fluctuations may change the proton spectrum only at energy substantially higher than $`E=4\times 10^{19}`$ eV. At this energy the half of energy losses is caused by $`e^+e^{}`$ production which does not fluctuate. Up to energy $`1\times 10^{20}`$ eV the photoproduction of pions occurs at the threshold in collisions with photons from high-energy tail of the Planck distribution, and fraction of energy lost does not fluctuate, being fixed by the threshold value. Indeed, for $`E_p=1\times 10^{20}`$ eV the minimal energy of CMB photon needed for pion production is $`ϵ=3\times 10^3`$ eV to be compared with energy of photon in the Planck distribution maximum $`ϵ_m=3.7\times 10^4`$ eV. The only fluctuating value is the interaction length.

The noticeable effect of fluctuations is expected for protons with energies $`E>1\times 10^{20}`$ eV.

As it is well known Lifshitz and Pitaevskii (1981); Ginzburg and Syrovatskii (1964), the kinetic equations give an adequate method to account for the fluctuations in interaction. Neglecting the conversion of proton to neutron (neutron decays back to proton with small energy loss) the kinetic equation for UHE protons with adiabatic energy losses and with $`p+\gamma p+e^++e^{}`$ and $`p+\gamma N+pions`$ scattering in collisions with CMB photons can be written down as follows:

$`{\displaystyle \frac{n_p(E,t)}{t}}=3H(t)n_p(E,t)+{\displaystyle \frac{}{E}}\left\{\left[H(t)E+b_{pair}(E,t)\right]n_p(E,t)\right\}`$ (16)
$`P(E,t)n_p(E,t)+{\displaystyle _E^{E_{max}}}𝑑E^{}P(E^{},E,t)n_p(E^{},t)+Q_{gen}(E,t),`$
where $`n_p(E,t)`$ is the number density of UHE protons per unit energy, $`Q_{gen}(E,t)`$ is the generation rate, given by Eq. (10) with $`m=0`$ and $`H(t)`$ is the Hubble parameter. The first term in the r.h.s. of Eq. (16) describes expansion of the universe. The energy loss $`b_{pair}(E)`$ due to $`e^+e^{}`$-pair production is treated as continuous energy loss. The photopion collisions are described with help of probability $`P(E,t)`$ for proton exit from energy interval $`(E,E+dE)`$ due to $`p\gamma `$-collisions and with help of their regeneration in the same energy interval described by probability $`P(E^{},E,t)`$. These two probabilities describe fluctuations in the interaction length and in fraction of energy lost in the interaction; the interaction length is equivalent in this picture to time of proton exit from energy interval $`dE`$.

The exit probability $`P(E,t)`$ due to collisions with CMB photons with temperature $`T`$ can be written as:

$`P(E,t)={\displaystyle \frac{cT}{4\pi ^2E^2}}{\displaystyle _{m_\pi +m_p}^{\mathrm{}}}𝑑E_cE_c(E_c^2m_p^2)\sigma (E_c)\mathrm{ln}\left[1\mathrm{exp}\left({\displaystyle \frac{E_c^2m_p^2}{4ET}}\right)\right],`$ (17)
where $`E_c`$ is the total c.m.s. energy of colliding proton and photon, $`E_c^{min}=m_\pi +m_p`$, $`\sigma (E_c)`$ is the photopion cross-section and $`T`$ is the CMB temperature at cosmological epoch $`t`$.

Similarly, in regeneration term of Eq. (16) the probability for a proton with energy $`E^{}`$ to produce a proton with energy $`E`$ is given by

$`P(E^{},E,t)={\displaystyle \frac{cT}{4\pi ^2(E^{})^2}}{\displaystyle \underset{E_c^{min}(x)}{\overset{\mathrm{}}{}}}𝑑E_cE_c(E_c^2m_p^2){\displaystyle \frac{d\sigma (E_c,E^{},E)}{dE}}\mathrm{ln}\left[1\mathrm{exp}\left({\displaystyle \frac{E_c^2m_p^2}{4E^{}T}}\right)\right],`$ (18)
where $`E^{}`$ and $`E`$ are the energies of primary and secondary protons, respectively, and $`x=E/E^{}`$. The minimum value of the allowed c.m.s. energy in this case is given by

$$E_c^{min}(x)=[m_\pi ^2/(1x)+m_p^2/x]^{1/2}.$$
(19)
This bound corresponds to the process with minimum invariant mass, namely to $`p+\gamma \pi ^0+p`$.

We have solved Eq. (16) numerically. The calculated spectrum $`J_p(E)=(c/4\pi )n_p(E,t_0)`$, where $`t_0`$ is the age of the universe at red-shift $`z=0`$, is presented in Fig. 15 as modification factor $`\eta (E)=J_p(E)/J_p^{\mathrm{unm}}(E)`$ for generation spectrum with $`\gamma _g=2.7`$ and $`E_{max}=1\times 10^{23}`$ eV (left panel) and $`E_{max}=1\times 10^{21}`$ eV (right panel). For comparison the modification factors for universal spectrum with continuous energy losses is also shown. The difference in these two spectra at highest energies must be due to fluctuations in energy losses, though formally we have to say that this is the difference between solution to kinetic equation (16) and the continuous energy loss approximation. For $`E_{max}=1\times 10^{23}`$ eV one can see the difference in the spectra about 25% at highest energies and tiny difference above intersection of $`\eta _{ee}`$ and $`\eta _{\mathrm{tot}}`$ curves. For $`E_{max}=1\times 10^{21}`$ eV the difference is small.

Note, that modification factors do not vanish at $`E_{\mathrm{max}}`$ even when generation function goes abruptly to zero, since both solutions vanish keeping the same value of ratios $`\eta (E)=J_p(E)/J_{\mathrm{unm}}(E)`$. It is easy to demonstrate analytically that ratio of flux in continuous loss approximation $`J_{\mathrm{cont}}(E)`$ to unmodified flux $`J_{\mathrm{unm}}(E)`$, given by Eq. (13) tends to $`H_0/\beta _{\mathrm{coll}}(E_{\mathrm{max}})`$, when $`EE_{\mathrm{max}}`$. From Fig. 15 one can see that this ratio coincides exactly with our numerical calculations (e.g. the analytical value is $`2.45\times 10^3`$ for $`E_{\mathrm{max}}=1\times 10^{23}`$ eV), and this gives a proof that our numerical calculations are correct.

The effect of interaction fluctuations is usually taken into account with help of Monte-Carlo simulations. The method of kinetic equation corresponds to averaging over large number of Monte-Carlo simulations, and if all other assumptions are the same, the results must coincide exactly. These assumptions include $`E_{\mathrm{max}}`$ and parameters of $`p\gamma `$-interaction. However, the existing Monte-Carlo simulations in most cases include some other assumptions in comparison with kinetic equations, which modify spectrum stronger than interaction fluctuations. One of them is discreteness in the source distribution (in kinetic equations the homogeneous distribution is assumed), the other is fluctuations of distances to the nearby sources.

It is possible however to make the comparison with Monte-Carlo simulations for homogeneously distributed sources and using identical interaction model. Such a comparison is discussed in Appendix C.

As to results presented here, it is necessary to emphasize that the difference of kinetic-equation solution and continuous-energy-loss approximation presented in Fig. 15 includes fluctuations, but not only fluctuations. The transition of kinetic equation to continuous-energy-loss equation depends on some other conditions which can fail. What the presented calculations demonstrate (see Fig. 15) is that continuous-energy-loss approximation describes with very good accuracy more reliable kinetic-equation solution, which in particular includes interaction fluctuations.

## VI  Transition from extragalactic to galactic cosmic rays

In the analysis above we obtained several indications that transition from extragalactic to galactic cosmic rays occurs at $`E1\times 10^{18}`$ eV. These evidences are summarized in Fig. 16.

The predicted spectrum above $`1\times 10^{18}`$ eV describes perfectly well the observed spectra: see the modification factor in the upper-left panel of Fig. 16 compared with AGASA data and in Fig. 8 with HiRes. However, at $`E<1\times 10^{18}`$ experimental modification factor becomes $`\eta >1`$, in contrast to definition $`\eta 1`$. It signals the appearance of a new component, which can be nothing but the galactic cosmic rays.

In the right panel the spectrum for rectilinear propagation from the sources with different separation $`d`$ and with $`\gamma _g=2.7`$ is compared with AGASA data. One can see that at $`E<1\times 10^{18}`$ eV the calculated extragalactic spectrum becomes less than that observed.

In the lower panel the similar comparison is shown for diffusive propagation with different diffusion regimes at lower energies Aloisio and Berezinsky (2005); Lemoine (2005). The random magnetic field with basic scale $`l_c=1`$ Mpc and magnetic field on this scale $`B_0=1`$ nG is assumed. The dash-dotted curve (universal spectrum) corresponds to the case when the separation between sources $`d0`$.

In all cases the transition from extragalactic to galactic component begins at $`E_b1\times 10^{18}`$ eV, with index $`b`$ for “beginning”.

What is the reason of this universality? We study the transition, moving from high towards low energies. $`E_b`$ is the beginning of transition (or its end, if one moves from low energies). $`E_b`$ is determined by energy $`E_{\mathrm{eq1}}=2.37\times 10^{18}`$ eV, where adiabatic and pair-production energy losses become equal. The quantitative analysis of this connection is given in Aloisio and Berezinsky (2005). We shall give here the semi-quantitative explanation.

The flattening of the spectrum occurs at energies $`EE_b`$, where $`E_b=E_{\mathrm{eq}}/(1+z_{\mathrm{eff}})^2`$ and $`z_{\mathrm{eff}}`$ should be estimated as red-shift up to which the main contribution to unmodified spectrum occurs. The simplified analytic estimate for $`\gamma _g=2.62.8`$ gives $`1+z_{\mathrm{eff}}1.5`$ and hence $`E_b1\times 10^{18}`$ eV. In fact, the right and lower panels of Fig. 16 present the exact calculations of this kind.

In experimental data the transition is searched for as a feature started at some low energy $`E_{2\mathrm{k}\mathrm{n}}`$ \- the second knee. Its determination depends on experimental procedure, and all we can predict is $`E_{2\mathrm{k}\mathrm{n}}<E_b`$. Determined in different experiments $`E_{2kn}(0.40.8)\times 10^{18}`$ eV.

The transition at the second knee appears also in the study of propagation of cosmic rays in the Galaxy (see e.g. P. L. Biermann *et al.* ; Hoerandel (2003); S. D. Wick, C. D. Dermer and Atoyan (2004)).

Being thought of as purely galactic feature, the position of the second knee in our analysis appears as direct consequence of extragalactic proton energy losses.

The transition at the second knee is illustrated by Fig. 17. The clue to understanding of this transition is given by the KASCADE data K.-H. Kampert *et al.* \[KASCADE collaboration\] (2001); J. R. Hoerandel *et al.* \[KASCADE collaboration\] (2002). They confirm the rigidity model, according to which position of a knee for nuclei with charge Z is connected with the position of the proton knee $`E_p`$ as $`E_Z=ZE_p`$. There are two versions of this model. One is the confinement-rigidity model (bending above the knee is due to insufficient confinement in galactic magnetic field), and the other is acceleration-rigidity model ($`E_{\mathrm{max}}`$ is determined by rigidity). In both models the heaviest nuclei (iron) start to disappear at $`E>E_{\mathrm{Fe}}=6.5\times 10^{16}`$ eV, if the proton knee is located at $`E_p2.5\times 10^{15}`$ eV. The shape of the spectrum above the iron knee ($`E>E_{\mathrm{Fe}}`$) is model dependent, with two reliably predicted features: it must be steeper than the spectrum below the iron knee ($`E<E_{\mathrm{Fe}}`$), shown by the dash-dot curve, and iron nuclei must be the dominant component there (see Fig. 17). The high energy part of the spectrum has a characteristic energy $`E_b`$, below which the spectrum becomes more flat, i.e. drops down when multiplied to $`E^{2.5}`$ (see Fig. 17). This part of the spectrum is shown for the diffusive propagation described in section VII.2. These two falling parts of the spectrum inevitably intersect at some energy $`E_{\mathrm{tr}}`$, which can be defined as transition energy from galactic to extragalactic cosmic rays. The ’end’ of galactic cosmic rays $`E_{\mathrm{Fe}}=6.5\times 10^{16}`$ eV and the beginning of full dominance of extragalactic component $`E_b1\times 10^{18}`$ eV differ by an order of magnitude. Note, that power-law extrapolation of the total galactic spectrum, shown by dot-dash line, beyond the iron knee $`E_{\mathrm{Fe}}`$ has no physical meaning in the rigidity models and must not be discussed.

The second-knee transition gives an alternative possibility in comparison with ankle-transition hypothesis known from end of 1970s. It is inspired by flattening of the spectrum at $`E_a1\times 10^{19}`$ eV seen in the AGASA and Yakutsk data (left panels in Fig. 8) and possibly at $`(0.51)\times 10^{19}`$ eV in the Hires data (the right panel in Fig. 8). Being multiplied to factor $`E^{2.5}`$, as in Fig. 17, the ankle transition looks very similar to that at the second knee. Note that in the latter case the ankle is just an intrinsic part of the dip.

The ankle transition has been recently discussed in Refs. Hillas (2004); Wibig and Wolfendale (2005); D. Allard, E. Parizot and Olinto (a, b); D. De Marco and Stanev (2005); Hillas .

In the ankle model it is assumed that galactic cosmic ray spectrum has a power-law shape $`E^\gamma `$ from the proton knee $`E_p2.5\times 10^{15}`$ eV to about $`E_a1\times 10^{19}`$ eV where it becomes steeper and crosses the more flat extragalactic spectrum (see the right panel of Fig. 17). The ankle transition in Fig. 17 is shown for extragalactic proton spectrum with generation index $`\gamma _g=2`$, while galactic spectrum, given by curve “gal.CR” is calculated as difference of the observed total spectrum and calculated extragalactic proton spectrum.

The ankle model has the problems with galactic component of cosmic rays. The spectrum at $`1\times 10^{18}1\times 10^{19}`$ eV is taken ad hoc to fit the observations, while in the second knee model this part of the spectrum is calculated with excellent agreement with the data. In the rigidity models the heaviest nuclei (iron) start to disappear at $`E>E_{\mathrm{Fe}}=6.5\times 10^{16}`$ eV. How the gap between $`1\times 10^{17}`$ eV and $`1\times 10^{19}`$ eV is filled?

Galactic protons start to disappear at $`E>2.5\times 10^{15}`$ eV. Where they came from at $`E>1\times 10^{17}`$ eV to be seen e.g. in the Akeno detector with fraction $`10\%`$?

The ankle model needs acceleration by galactic sources up to $`1\times 10^{19}`$ eV (at least for iron nuclei), which is difficult to afford. The second knee model ameliorates this requirement by one order of magnitude.

The second-knee model predicts the spectrum shape down to $`1\times 10^{18}`$ eV with extremely good accuracy ($`\chi ^2`$/d.o.f.= 1.12 for Akeno-AGASA and $`\chi ^2`$/d.o.f.= 1.03 for HiRes). In the ankle model one has to consider this agreement as accidental, though such hypothesis has very low probability, determined by $`\chi ^2`$ cited above. As an alternative the ankle model-builders can suggest only hopes for future development of galactic propagation models to be as precisely calculated as the dip.

## VII  Astrophysical sources of UHECR

In the sections above we have performed the model-independent analysis of spectra of extragalactic protons interacting with CMB. We have calculated the features of the proton spectrum assuming the power-law generation spectrum $`E^{\gamma _g}`$ valid at $`E1\times 10^{18}`$ eV, and compared predicted features with observations. We found that proton dip, a model-independent feature at energy between $`1\times 10^{18}`$ eV and $`4\times 10^{19}`$ eV, is well confirmed by observations. Only two free parameters are involved in fitting of observational data: $`\gamma _g=2.7`$ (the allowed range is 2.55 - 2.75) and the flux normalization constant. The various physical phenomena included in calculations, such as discreteness in the source distribution, the different modes of propagation (rectilinear and diffusive), cosmological evolution with parameters similar to AGN evolution, fluctuations in $`p\gamma `$ interaction etc, do not upset this agreement.

The transition of extragalactic to galactic cosmic rays is also discussed basically in model-independent manner.

In this Section we shall discuss the models: realistic energy spectra, the sources and the models for transition from extragalactic to galactic cosmic rays.

The UHECR sources have to satisfy two conditions: they must be very powerful and must accelerate particles to large $`E_{\mathrm{max}}>1\times 10^{21}`$ eV. There is one more restriction S. L. Dubovsky, P. G. Tinyakov and Tkachev (2000); Fodor and Katz (2001); H. Yoshiguchi, S. Nagataki and Sato (2003); G. Sigl, F. Miniati and Ensslin (2003a); H. Yoshiguchi, S. Nagataki and Sato (2004a); P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a), coming from observation of small-scale clustering: the space density of the sources should be $`(13)\times 10^5`$ Mpc-3 probably with noticeable uncertainty in this value. Thus, these sources are more rare, than typical representatives of AGN, e.g. the Seyfert galaxies, whose space density is $`3\times 10^4`$ Mpc-3. The sources could be the rare types of AGN, and indeed the analysis of Tinyakov and Tkachev (2001) show statistically significant correlation between directions of particles with energies $`(48)\times 10^{19}`$ eV and directions to AGN of the particular type – BL Lacs (see also criticism N. W. Evans, F. Ferrer and Sarkar (2004) and reply Tinyakov and Tkachev (2004)). The acceleration in AGN can provide the maximum energy of acceleration up to $`10^{21}`$ eV for non-relativistic shock acceleration (see e.g. Biermann and Strittmatter (1987)).

The relativistic shock acceleration can occur in AGN jets. Acceleration to $`E_{\mathrm{max}}10^{21}`$ eV in the AGN relativistic shocks is questionable (see discussion below).

Gamma Ray Bursts (GRBs) are another potentially possible sources of UHECR Milgrom and Usov (1995); Vietri (1995); Waxman (1995). They have very large energy output and can accelerate particles up to $`10^{21}`$ eV Waxman (1995); Vietri (1995). These sources have, however, the problems with explaining small-angle anisotropy and with energetics (see discussion below).

### VII.1  Spectra

The assumption of the power-law generation spectrum with $`\gamma _g=2.7`$ extrapolated to $`E_{\mathrm{min}}1`$ GeV results in too large emissivity required for observed fluxes of UHECR. To avoid this problem the broken generation spectrum has been suggested in Refs. V. Berezinsky, A. Z. Gazizov and Grigorieva (a, b):

$$Q_{\mathrm{gen}}(E)=\{\begin{array}{cc}E^2\hfill & \mathrm{at}EE_c\hfill \\ E^{2.7}\hfill & \mathrm{at}EE_c\hfill \end{array},$$
(20)
where $`Q_{\mathrm{gen}}(E)`$ is the generation function (rate of particle production per unit of comoving volume), defined by Eqs. (9) and (10), and $`E_c`$ was considered as a free parameter.

Recently it was demonstrated that broken generation spectrum can naturally emerge under most reasonable physical assumption. In Ref. Kachelriess and Semikoz (2005b) it has been argued that while spectrum $`E^2`$ is universal for non-relativistic shock acceleration, the maximum acceleration energy $`E_{\mathrm{max}}`$ is not, being dependent on the physical characteristics of a source, such as its size, regular magnetic field etc. Distribution of sources over $`E_{\mathrm{max}}`$ results in steepening of the generation function, so that the distribution of the sources $`dn/dE_{\mathrm{max}}E_{\mathrm{max}}^\beta `$ explains the observational data, if $`\beta =1.51.6`$ Kachelriess and Semikoz (2005b).

In Appendix E we address the generalized problem: what should be the distribution of spectral emissivity over $`E_{\mathrm{max}}`$ to provide the generation function with the broken spectrum. We use there the notation $`\epsilon E_{\mathrm{max}}`$ and introduce the spectral emissivity $`(\epsilon )=n_s(\epsilon )L_p(\epsilon )`$, where $`L_p(\epsilon )`$ is particle luminosity of a source and $`n_s(\epsilon )n_s(\epsilon ,L_p(\epsilon ))`$ is the space density of the sources. The total emissivity is given by $`_0=(\epsilon )𝑑\epsilon `$. Distribution of spectral emissivity $`(\epsilon )`$ over maximal energies $`\epsilon `$ determines the energy steepening of the generation function $`Q_{\mathrm{gen}}(E)`$ at energy $`E=\epsilon _{\mathrm{min}}`$ in the distribution. This function is calculated analytically for arbitrary $`(\epsilon )`$, assuming that $`\epsilon `$ is confined to the interval ($`\epsilon _{\mathrm{min}},\epsilon _{\mathrm{max}})`$. It is demonstrated that $`Q_{\mathrm{gen}}(E)`$ can be the power-law function $`E^{\gamma _g}`$ exactly, only if $`(\epsilon )`$ distribution is power-law, too $`(\epsilon ^\beta )`$. For the source generation function $`q_{\mathrm{gen}}(E)E^{(2+\alpha )}`$ at $`E\epsilon `$, the generation index in the interval $`\epsilon _{\mathrm{min}}E\epsilon _{\mathrm{max}}`$ is found to be $`\gamma _g=1+\alpha +\beta `$, including the case $`\alpha =0`$. The steepening of the generation spectrum from $`2+\alpha `$ to $`\gamma _g`$ occurs approximately at energy $`E_c=\epsilon _{\mathrm{min}}`$. At energy $`E>\epsilon _{\mathrm{max}}`$ the spectrum is suppressed as $`\mathrm{exp}(E/\epsilon _{\mathrm{max}})`$ or stronger.

In the applications we are interested in two cases: $`\alpha =0`$ (non-relativistic shocks) and $`\alpha =0.20.3`$ (ultra-relativistic shocks). In the latter case the term $`E^2`$ in Eq. (20) should be substituted by $`E^{(2+\alpha )}`$. The energy $`E_c=\epsilon _{min}`$ in Eq. (20) is considered as a free parameter.

### VII.2 Active Galactic Nuclei

The AGN as sources of UHECR meet the necessary requirements: (i) to accelerate particles to $`E_{\mathrm{max}}10^{21}`$ eV, (ii) to provide the necessary energy output and (iii) to have the space density $`n_s(13)\times 10^5`$ Mpc-3, required by small-scale clustering. We shall discuss below these problems in some details.

#### VII.2.1 Acceleration and spectra

The flow of the gas in AGN jet can be terminated by the non-relativistic shock which accelerates protons or nuclei in the radio lobe up to $`E_{\mathrm{max}}10^{21}`$ eV with spectrum $`E^2`$ Biermann and Strittmatter (1987).

In some cases the observed velocities in AGN jets are ultra-relativistic with Lorentz factor up to $`\mathrm{\Gamma }510`$. It is natural to assume there the existence of internal and external ultra-relativistic shocks. Acceleration in relativistic shocks relevant for UHECR has been recently studied in Refs. Vietri (1995); Waxman (1995); Gallant and Achterberg (1999); J. G. Kirk, A. W. Guthmann, Y. A. Gallant and Achterberg (2000); Vietri (2003); Blasi and Vietri (2005); Lemoine and Pelletier (2003); Lemoine and Revenu (2006). The acceleration spectrum is $`E^{\gamma _g}`$ and in the case of isotropic scattering of particles upstream and downstream, the spectrum index is $`\gamma _g=2.23\pm 0.01`$ J. G. Kirk, A. W. Guthmann, Y. A. Gallant and Achterberg (2000). However, recently it was understood that this result depends on scattering properties of the medium Blasi and Vietri (2005); Lemoine and Revenu (2006), and the spectrum can be steeper. In the regime of large angle scattering in Blasi and Vietri (2005), $`\gamma _g=2.7`$ was found possible for the shock with velocity $`u(0.80.9)c`$ and compression ratio $`r=2`$. In Monte-Carlo simulation Lemoine and Revenu (2006) it is demonstrated that effect of compression of upstream magnetic field results in increasing of $`\gamma _g`$ up to the limiting value 2.7 in ultra-relativistic case $`\mathrm{\Gamma }_{\mathrm{sh}}1`$.

The maximal energy of acceleration $`E_{\mathrm{max}}`$ is a controversial issue. While in most works (very notably Gallant and Achterberg (1999)) it is obtained that $`E_{\mathrm{max}}`$ cannot reach $`10^{21}`$ eV for all realistic cases of relativistic shocks, the authors of Ref. M. Vietri, D. De Marco and Guetta (2003) argue against this conclusion.

We shall divide a problem with $`E_{\mathrm{max}}`$ into two: the reliably estimated energy gain in relativistic shocks and model-dependent absolute value of $`E_{\mathrm{max}}`$. The energy gain in the first full cycle of particle reflection upstream-downstream-upstream ($`udu`$) is about $`\mathrm{\Gamma }_{\mathrm{sh}}^2`$. The next reflections are much less effective, as was first observed in Gallant and Achterberg (1999): a particle lives a short time in the upstream region, before it is caught up by the shock. As a result, a particle is deflected by upstream magnetic field only to a small angle, and thus it occurs in the downstream region with approximately the same energy as in the first cycle. Then the energy upstream will be also almost the same as in the first cycle. According to Monte-Carlo simulation Lemoine and Revenu (2006) the average energy gain per each successive $`udu`$ cycle is only 1.7 (see Fig. 4 in Lemoine and Revenu (2006) with clear explanation). With these energy gains ($`\mathrm{\Gamma }_{\mathrm{sh}}^2`$ in the first $`udu`$ cycle and with $`2`$ in each successive cycle) it is not possible to get $`E_{\mathrm{max}}10^{21}`$ eV in the conservative approach for AGN and GRBs, but there are some caveats in this conclusion as indicated in Gallant and Achterberg (1999); M. Vietri, D. De Marco and Guetta (2003). The magnetic field in the upstream region can be large, and then the deflection angle of a particle after the shock crossing is large, too. Relativistic shock acceleration can operate in medium filled by pre-accelerated particles, and thus initial energy can be high.

There are some other mechanisms of acceleration to energies up to $`10^{21}`$ eV relevant for AGN: unipolar induction and acceleration in strong electromagnetic waves (see V. S. Berezinsky, S. V. Bulanov, V. A. Dogiel, V. L. Ginzburg and Ptuskin (1990) for description and references). The mechanisms of jet acceleration have very special status.

The observed correlations between arrival directions of particles with energies $`(48)\times 10^{19}`$ eV and BL Lacs Tinyakov and Tkachev (2001) imply the jet acceleration. This is because BL Lacs are AGN with jets directed towards us. For this correlation the propagation of particles (most probably protons) with energies above $`4\times 10^{19}`$ eV must be (quasi)rectilinear, that can be realized in magnetic fields found in MHD simulations in K. Dolag, D. Grasso, V. Springel and Tkachev (2004) (see however the simulations in G. Sigl, F. Miniati and Ensslin (2003b); G. Sigl, F. Miniati and Ensslin (2004) with quite different results).

An interesting mechanism of jet acceleration, called pinch acceleration, was suggested and developed in plasma physics V. V. Vlasov, S. K. Zhdanov and Trubnikov (1990). It is based on pinch instability well known in plasma physics, both theoretically and observationally. The electric current along jet produces the toroidal magnetic field which stabilizes the jet flow. The pinch instability is caused by squeezing the tube of flow by magnetic field of the current. It results in increasing the electric current density and magnetic field. The magnetic field compresses further the tube and thus instability develops. The acceleration of particles in V. V. Vlasov, S. K. Zhdanov and Trubnikov (1990) is caused by hydrodynamical increase of velocity of the flow and by longitudinal electric field produced in the pinch. This process is illustrated by Fig. 18.

The pinch acceleration has been developed for tokamaks and was confirmed there by observations. The generation spectrum is uniquely predicted as

$$q_{\mathrm{gen}}(E)E^{\gamma _g},\mathrm{with}\gamma _g=1+\sqrt{3}=2.73.$$
(21)
Scaling a size of the laboratory tube to the cosmic jet, one obtains $`E_{\mathrm{max}}`$ exceeding $`10^{22}`$ eV. The particle beam undergoes a few pinch occurrences during traveling along cosmic jet (a few kpc in case of AGN), and thus spectrum obtains a low energy cut at high value of $`E_{\mathrm{min}}`$. This mechanism needs more careful study.

We finalize this subsection concluding, that there are at least three possibilities for the broken generation spectrum with $`q_{\mathrm{gen}}E^{2.7}`$ at $`E>E_c`$: The acceleration by non-relativistic shocks with distribution of sources (more precisely emissivity) over $`E_{\mathrm{max}}`$, the acceleration by relativistic shocks, where $`E_c=E_{\mathrm{min}}`$ naturally occurs due to the first $`udu`$ cycle as $`\mathrm{\Gamma }_{\mathrm{sh}}^2E_{\mathrm{in}}`$ with $`E_{\mathrm{in}}`$ being particle energy before acceleration, and in the pinch acceleration where $`\gamma _g=1+\sqrt{3}`$ is rigorously predicted and $`E_{\mathrm{min}}`$ appears due to several pinch occurrences. The latter mechanism provides acceleration along a jet, necessary for correlations with BL Lacs.

#### VII.2.2 Spectra of UHECR from AGN

We discuss here the diffuse energy spectra from AGN. They are model-dependent because one should specify the distance between sources $`d`$, the mode of propagation (e.g. rectilinear or diffusive) and the critical energy $`E_c`$ in the broken spectrum of generation. The universal spectrum is valid when characteristic distance $`d`$ between sources is smaller than other propagation lengths, most notably energy attenuation length $`l_{\mathrm{att}}(E)`$ and diffusion length $`l_{\mathrm{diff}}(E)`$. As was discussed in Section IV the dip, seen at energies $`1\times 10^{18}<E<4\times 10^{19}`$ eV, is the robust feature, which is not modified by any known phenomenon except presence of extragalactic nuclei with fraction higher than $`20\%`$. This fraction depends on model of acceleration and can be small for some acceleration mechanisms, e.g. for acceleration by relativistic shocks.

The prediction of the GZK cutoff is model dependent: the shape of this feature can be modified by different values of $`E_{\mathrm{max}}`$, by local overdensity/deficit of the sources (see Section III.1) and by discreteness in the source distribution (see Fig. 9). It can also be mimicked by the shape of acceleration cutoff near $`E_{\mathrm{max}}`$.

The spectrum predicted for $`E<1\times 10^{18}`$ eV is also model dependent. At these energies for all distances $`d`$ reasonable for AGN the universal spectrum is modified due to propagation in magnetic fields, because $`l_{\mathrm{diff}}`$ is small. In Fig. 17 (left panel) the spectrum for diffusion in random magnetic field is shown. At energies $`E>1\times 10^{18}`$ eV the spectrum remains universal.

We have to introduce in calculations the broken generation spectrum (20) to provide the reasonable luminosities of AGN. Normalized by the AGASA flux, the emissivity $`_0`$ is given by $`3.5\times 10^{46}`$ ergs Mpc-3yr-1, $`1.6\times 10^{47}`$ ergs Mpc-3yr-1 and $`7.0\times 10^{47}`$ ergs Mpc-3yr-1 for $`E_c`$ equals to $`1\times 10^{18}`$ eV, $`1\times 10^{17}`$ eV and $`1\times 10^{16}`$ eV, respectively.

To calculate the luminosities of AGN, $`L_p=_0/n_s`$, we take the density of the sources from small-scale clustering as $`n_s=2\times 10^5`$ Mpc-3 P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a). The luminosities are found to be quite reasonable for AGN: $`5.6\times 10^{43}`$$`2.5\times 10^{44}`$ and $`1.1\times 10^{45}`$ erg s-1 for $`E_c`$ equals to $`1\times 10^{18}`$ eV, $`1\times 10^{17}`$ eV and $`1\times 10^{16}`$ eV, respectively.

The calculated spectra are shown in Fig. 19 for propagation in magnetic field, as described in Fig. 10 for the Bohm diffusion in the low-energy regime, for $`d=50`$ Mpc, $`\gamma _g=2.7`$, $`E_c=1\times 10^{18}`$ eV and $`m=0`$, inspired by BL Lacs, which show no positive evolution. The spectra for AGN evolutionary models are shown in Fig. 14 in comparison with case $`m=0`$. For calculations of spectra with different values of $`E_c`$ see V. Berezinsky, S. I. Grigorieva and Hnatyk (2004).

#### VII.2.3 Transition to galactic cosmic rays

For the diffusive propagation the transition is shown in Fig. 17. In our calculations we followed Ref. Aloisio and Berezinsky (2005), considering the same case of diffusive propagation in random magnetic fields with basic coherent scale $`l_c=1`$ Mpc and $`B_0=1`$ nG on this scale. These parameters define the diffusion coefficient at very high energies, while at energies $`E<1\times 10^{18}`$ eV at interest we assume the Bohm diffusion. The distance between sources is fixed as $`d=50`$ Mpc. As Fig. 17 shows the iron galactic spectrum given by curve “gal.Fe” describes well the KASCADE flux at $`E1\times 10^{17}`$ eV and fits the Hall (drift) diffusion spectrum at $`E>1\times 10^{17}`$ eV. The predicted fraction of iron in the total flux is given in Table 1 as function of energy. This fraction is higher than in Aloisio and Berezinsky (2005) mostly because we used the HiRes data, instead of AGASA in Aloisio and Berezinsky (2005).

The detailed study of transition for rectilinear propagation for different values of $`E_c`$ has been performed in V. Berezinsky, S. I. Grigorieva and Hnatyk (2004). The fraction of iron-nuclei varies with $`E_c`$.

#### VII.2.4 Signatures of AGN

AGN are one of the best candidates for UHECR sources. They have high luminosities to provide the observed UHECR flux and have quite efficient mechanisms of acceleration. The predicted spectra are in good agreement with observations and transition to galactic cosmic rays is well described using diffusive or quasi-rectilinear propagation and broken generation spectrum. The density of sources $`n_s(13)\times 10^5`$ Mpc-3, found from small-scale clustering, corresponds well to space density of powerful AGN. However, in the diffuse spectrum we cannot predict any specific signature of AGN, and probably only observations of the direct flux of primary (e.g. protons) or secondary (e.g. neutrinos or gamma-rays) radiation can give the direct evidence for AGN as sources of UHECR.

In this respect the correlations with BL Lacs Tinyakov and Tkachev (2001, 2004) is the direct indication to AGN as UHECR sources. The protons can be the signal carriers. The important point is that in structured (perforated) magnetic field, there are directions (holes) with weak magnetic fields in which BL Lacs can be seen V. Berezinsky, A. Z. Gazizov and Grigorieva (b). In MHD simulation K. Dolag, D. Grasso, V. Springel and Tkachev (2004) it was found indeed that in many large scale structures like voids and filaments the magnetic field is weak and protons can propagate there with small deflections.

Correlation with BL Lacs imply the jet acceleration, because BL Lacs are AGN with jets directed to an observer. The pinch acceleration mechanism fits well the picture above, providing $`\gamma _g=2.7`$, high $`E_{\mathrm{max}}`$ and high $`E_{\mathrm{min}}`$, needed for reasonable luminosity $`L_p`$.

### VII.3 Gamma Ray Bursts

GRBs are another candidates for UHECR sources. Like AGN they have high energy output and can accelerate particles to UHEs, as first have been proposed in Milgrom and Usov (1995); Vietri (1995); Waxman (1995)

GRB starts with instant energy release $`W_{\mathrm{tot}}`$, which occurs most probably due to SN explosion. This energy is injected into a small volume with initial radius $`R_i10^7`$ cm in a form of a jet with small opening angle $`\theta `$. The extremely large energy release within a small volume results in very high temperature of the gas, which consists mostly of photons and electron-positron pairs with a small admixture of initially injected baryons. The pressure of relativistic $`e^+e^{}\gamma `$ gas accelerates fireball, its temperature falls down and the Lorentz factor rises. At later stages a fireball becomes the baryonic jet moving with Lorentz factor $`\mathrm{\Gamma }_0=W_{\mathrm{tot}}/M_bc^2100`$, where $`M_b`$ is initial baryon loading.

The jet is characterized according to observations D. A. Frail *et al.* (2001) by small opening angle $`\theta 4^{}`$ and by corresponding solid angle $`\mathrm{\Omega }=\pi \theta ^2`$. The external observer sees only part of the external fireball surface within the angle $`1/\mathrm{\Gamma }_0<\theta `$. Because if this fact, the convenient formalism is to consider the fictitious spherically-symmetric fireball with equivalent isotropic energy output $`W_{\mathrm{tot}}`$ and equivalent isotropic GRB rate $`\nu _{\mathrm{GRB}}`$ (in units of Mpc-3 yr-1). The real energy output and GRB rate are $`\stackrel{~}{W}_{\mathrm{tot}}=(\mathrm{\Omega }/4\pi )W_{\mathrm{tot}}`$ and $`\stackrel{~}{\nu }_{\mathrm{GRB}}=(4\pi /\mathrm{\Omega })\nu _{\mathrm{GRB}}`$, with the same emissivity $`=W_{\mathrm{tot}}\nu _{\mathrm{GRB}}`$.

According to observations the equivalent isotropic energy output $`W_{\mathrm{GRB}}1\times 10^{53}`$ erg, and spherically-symmetric GRB rate is $`\nu _{\mathrm{GRB}}0.5\times 10^9`$ Mpc-3yr-1 with emissivity $`_{\mathrm{GRB}}0.5\times 10^{44}`$ erg Mpc-3yr-1. The more detailed estimates differ considerably from this value: $`_{\mathrm{GRB}}0.6\times 10^{43}`$ erg Mpc-3yr-1 (M. Schmidt Schmidt (1999, 2001)) and $`_{\mathrm{GRB}}1.3\times 10^{44}`$ erg Mpc-3yr-1 (Frail *et al.* D. A. Frail *et al.* (2001)). Note, that difference in spectral intervals used in Schmidt (1999, 2001) and D. A. Frail *et al.* (2001) cannot explain alone the differences between the above-cited emissivities. We conclude thus that $`_{\mathrm{GRB}}(0.613)\times 10^{43}`$ erg Mpc-3yr-1, to be compared with higher estimate by Waxman Waxman (2001) $`_{\mathrm{GRB}}(330)\times 10^{43}`$ erg Mpc-3yr-1.

There are two sites of acceleration in a jet: mildly relativistic internal shocks inside a fireball Waxman (1995) and ultra-relativistic external shock Vietri (1995).

#### VII.3.1 Internal shock acceleration and UHECR

We shall describe here shortly the acceleration inside fireball, following Waxman (1995, 2001); Piran (1999).

Gamma-radiation is produced as synchrotron emission of the shock-accelerated electrons at the stage when fireball moves with constant Lorenz-factor $`\mathrm{\Gamma }_0=W_{\mathrm{tot}}/M_bc^2`$. Acceleration occurs due to internal shocks produced by motion of fireball sub-shells with different Lorentz factors $`\mathrm{\Gamma }`$. This phenomenon is assumed to be responsible for the observed GRBs and for their time substructure. In the observer frame the jet looks like a narrow shell with width $`\mathrm{\Delta }R/\mathrm{\Gamma }_0^2`$, where $`R`$ is a radius, i.e. distance from outer surface to the source (the GRB engine). It is assumed that this shell consists of mildly relativistic sub-shells, moving relatively each other in the rest frame of the fireball. The inner sub-shells move faster than the outer ones, and their collisions produce mildly relativistic shocks. The GRB spikes arise due to shell substructure. The range of distances $`R_\gamma `$ where GRB is produced can be estimated from the observed minimum duration of spikes, $`\tau _{\mathrm{sp}}>1`$ ms, and duration of GRB, $`\tau _{\mathrm{GRB}}110`$ s. Time between collisions of two sub-shells can be estimated in the observer frame Piran (1999) as $`\delta t\delta /v\mathrm{\Gamma }_0^2\delta /c`$, where $`\delta `$ is the distance between two colliding subshells. Thus, $`R_\gamma =c\delta t`$ is limited between $`R_\gamma ^{\mathrm{min}}\mathrm{\Gamma }_0^2c\tau _{\mathrm{sp}}3\times 10^{11}`$ cm and $`R_\gamma ^{\mathrm{max}}\mathrm{\Gamma }_0^2c\tau _{\mathrm{GRB}}3\times 10^{14}3\times 10^{15}`$ cm. The energetics is described by $`W_{\mathrm{GRB}}=f_{\mathrm{GRB}}W_{\mathrm{tot}}`$, with $`f_{\mathrm{GRB}}<0.5`$ at least, because half of the total energy is released in afterglow. More realistically $`f_{\mathrm{GRB}}0.2`$.

The exit of all particles from the jet occurs only through the front spherical surface, provided by condition $`\theta >1/\mathrm{\Gamma }_0`$.

Acceleration of electrons by internal shocks implies acceleration of protons. According to Waxman (1995, 2001) the generation spectrum $`E^2`$, $`E_{\mathrm{max}}1\times 10^{21}`$ eV, $`E_{\mathrm{min}}\mathrm{\Gamma }_0m_pc^2`$ and $`W_pW_{\mathrm{GRB}}`$ (in the observer frame).

As a matter of fact one must write the expressions above as $`E_{\mathrm{min}}f_a\mathrm{\Gamma }_0m_pc^2`$ and $`W_pf_aW_{\mathrm{GRB}}`$, where $`f_a`$ is a factor of adiabatic cooling. This phenomenon has been studied in Rachen and Meszaros (1998) and adiabatic energy losses were found to be very severe. As was indicated there the efficient mechanism, which diminishes the influence of adiabatic cooling, is the neutron mechanism of UHECR exit from expanding fireball. However, some restrictions obtained for this mechanism make its application rather limited. Waxman Waxman (2001) fought back this criticism by argument that production of HE protons in the shell at later periods corresponding to a distance $`R_\gamma ^{\mathrm{max}}`$ eliminates the problem of adiabatic cooling. We shall argue below that this problem exists.

There are three stages of adiabatic energy losses. The first one occurs during acceleration process and it operates at distances $`R_\gamma ^{\mathrm{min}}<R<R_\gamma ^{\mathrm{max}}`$. The adiabatic energy losses are described in the standard shock-acceleration equation by the term $`(p/3)div\stackrel{}{u}f(r,t,p)/p`$, where $`\stackrel{}{u}`$ is hydrodynamical velocity of the gas flow and $`f(r,t,p)`$ is the distribution function. In our case the adiabatic energy losses are large due to expansion of the shell. Adiabatic energy losses diminish the fraction of the shock energy transferred to accelerated particles, which results in relation $`W_p=f_a^{(1)}W_{\mathrm{GRB}}`$, where $`f_a^{(1)}<1`$ is adiabatic factor. Qualitatively it can be explained in the following way. The first portion of accelerated particles would suffer severe adiabatic energy losses $`f_aR_\gamma ^{\mathrm{min}}/R_\gamma ^{\mathrm{max}}`$, if further acceleration were absent. In fact, these particles are re-accelerating and loosing energy simultaneously, as corresponding diffusion equation describes it. This process diminishes only the fraction of energy $`f_a`$ transferred to accelerated particles. It can be estimated as $`f_aR_\gamma ^{\mathrm{eff}}/R_\gamma ^{\mathrm{max}}`$, but for accurate estimate the solution of diffusion equation is needed. Unfortunately, this model exists only in the form of estimates.

The second stage starts after acceleration ceases. The end of GRB radiation signals this moment. Magnetic field in the expanding shell slowly diminishes, accelerated electrons radiate away their energy fast, but protons are still confined in the shell, forming relativistic gas and loosing energy adiabatically. Duration of this stage is different for protons of different energies and the spectrum $`E^2`$ is distorted. The protons with $`10^{19}10^{21}`$ eV are the first to accomplish this stage and enter the regime of free expansion (stage 3). This energy range is most interesting for UHECR . The adiabatic factor $`f_a^{(2)}`$ depends on the time of diminishing of equipartition magnetic field. The detailed calculations are needed to evaluate $`f_a^{(2)}`$ for the range $`10^{19}10^{21}`$ eV.

At the third stage when the shell of the protons with energies $`10^{19}10^{21}`$ eV propagates quasi-rectilinearly in the external magnetic field, the energy losses are also present. They are caused by the collective effect of transferring proton momenta to the macroscopic volume of gas due to scattering in magnetic field.

As a toy model, let us consider the spherically symmetric shell of protons with energies $`10^{19}10^{21}`$ eV, and make an order of magnitude estimates. At the distance of order of gyroradius $`r_L`$, a proton transfers half of its initial (radial) momentum to the gas. Since protons propagate in the radial direction almost with speed of light, it follows from causality that volume of gas, which absorbs the lost momentum cannot be larger than $`r_L^3`$. The equation of radial momentum balance reads

$$0.5W_p/c\mathrm{\Gamma }r_L^3\rho _gc,$$
(22)
where $`\rho _g10^{24}`$ g/cm3 is the density of the ambient gas and $`\mathrm{\Gamma }`$ is Lorentz factor, obtained by this gas volume. For the spectrum $`E^2`$ the total energy of protons with $`10^{19}E10^{21}`$ eV is $`W_p0.2W_{\mathrm{GRB}}`$. For magnetic field in the presupernova wind at distance $`r10^{16}`$ cm we accept the value of magnetic field $`B_01`$ G, following the estimate of M. Vietri, D. De Marco and Guetta (2003) $`0.110`$ G (see also Waxman (2001)). Then for $`W_{\mathrm{GRB}}1\times 10^{53}`$ erg we obtain $`\mathrm{\Gamma }10^5`$ and for $`W_{\mathrm{GRB}}5\times 10^{50}`$ erg (jet value) $`\mathrm{\Gamma }10^3`$. In fact for the jet case the effective mass of the target is much less and $`\mathrm{\Gamma }`$ is larger. Such gas shell produces the relativistic shock and considered scenario radically changes.

If considered UHE protons capture electrons from the media with its frozen magnetic field, the pressure of the beam to the gas further increases.

In fact, we considered above only one mechanism for energy transfer to the gas. It is clear, that nature of this phenomenon is very general: the energy density of proton beam is so large, that part of its energy is transferred to magnetized ambient plasma. Only when this density drops to some critical value and collective effects disappear, the protons can propagate as individual particles, not exciting the ambient plasma.

Therefore, the total adiabatic factor $`f_a`$ cannot be of order of 1.

In the left panel of Fig. 20 the calculated flux of UHE protons accelerated by internal shocks is presented by curve 1. It corresponds to generation spectrum $`E^2`$ with $`E_{\mathrm{max}}=1\times 10^{21}`$ eV and $`E_{\mathrm{min}}=\mathrm{\Gamma }_0m_pc^2=1\times 10^{11}`$ eV. This is the case of absence of adiabatic energy losses $`f_a=1`$. The spectrum is calculated in the standard way taking into account the energy losses on CMB radiation. The emissivity needed to fit the observational data is $`_{\mathrm{CR}}=2.2\times 10^{45}`$ erg Mpc-3yr-1. It must be equal to $`W_p\nu _{\mathrm{GRB}}`$ and thus it should not exceed $`_{\mathrm{GRB}}`$ given above. As a matter of fact it is by factor 20 - 300 higher. If adiabatic energy losses for protons with energies $`10^{19}10^{21}`$ eV is $`f_a(>1\times 10^{19})`$, the predicted flux in this energy range is additionally suppressed by this factor and discrepancy increases correspondingly.

If to assume $`f_a(E1\times 10^{19}\mathrm{eV})1`$, while $`f_a(E1\times 10^{19}\mathrm{eV})1`$, the energy balance is not changed, because the total energy (including $`W_p`$ and that adiabatically lost) remains the same.

In the right panel of Fig. 20 we show the spectrum with the assumed distribution of sources over $`E_{\mathrm{max}}`$ like in case of AGN. In this case the required emissivity is higher $`_{\mathrm{CR}}=3.5\times 10^{46}`$ erg Mpc-3yr-1.

These results have been already presented in V. Berezinsky, A. Z. Gazizov and Grigorieva (a). Waxman Waxman (2001) did not understand that we considered two different cases (now being put separately in two different panels in Fig. 20) and argued in fact only against the case shown in the right panel of Fig. 20. The left panel shows the same case as considered by Waxman. Our discrepancy in emissivity is explained as follows. We normalize our calculations by emissivity $`_p`$ (the energy released in cosmic ray protons per unit comoving volume and per unit time) and compare it with $`_{\mathrm{GRB}}=W_{\mathrm{GRB}}\nu _{\mathrm{GRB}}`$. Waxman normalizes his flux by value $`E^2\dot{n}_{CR}/dE=0.8\times 10^{44}`$ erg/Mpc3yr Waxman (2001). To obtain emissivity, this value should be multiplied by $`\mathrm{ln}(E_{\mathrm{max}}/E_{\mathrm{min}})`$, which is equal to 23 for $`E_{\mathrm{max}}=1\times 10^{21}`$ eV and $`E_{\mathrm{min}}=1\times 10^{11}`$ eV. The obtained UHECR emissivity roughly agrees with our value, but it exceeds by an order of magnitude what Waxman cites in Waxman (2001) as the highest value of GRB emissivity ($`3\times 10^{44}`$ erg Mpc-3yr-1). Note that the comparison is made for unrealistic case of absence of adiabatic energy losses.

#### VII.3.2 External shock acceleration and UHECR

The flow of ultra-relativistic baryonic jet is terminated by the external ultra-relativistic shock propagating in pre-supernova wind. The acceleration by such shocks has been described (with relevant references) in subsection VII.2.1. As a shock propagates in interstellar medium, the protons, experiencing the different number of $`udu`$ cycles, move together with the shock ahead of it (upstream) or behind it (downstream). The “fresh” particles, which undergo one or a few ($`n`$) cycles have energy $`E\mathrm{\Gamma }_{\mathrm{sh}}^22^nE_i`$, the “old” particles approach $`E_{\mathrm{max}}eBR_{\mathrm{sh}}\mathrm{\Gamma }_{\mathrm{sh}}`$ Gallant and Achterberg (1999), where $`E_i`$ is initial proton energy upstream, $`R_{\mathrm{sh}}`$ is a radius of the shock (the distance from GRB origin), and $`\mathrm{\Gamma }_{\mathrm{sh}}100`$ is the Lorentz-factor of the shock. The acceleration time $`t_aE/(ecB\mathrm{\Gamma }_{\mathrm{sh}})`$ Gallant and Achterberg (1999) is smaller than the age of the shock $`tR_{\mathrm{sh}}/c`$ for all particles except those with $`EE_{\mathrm{max}}`$.

At this stage of our study we are interested only in the total energy transferred to accelerated protons, disregarding what happens to these particles later. In other words we are interested only in energy balance, using the common assumption $`W_pW_{\mathrm{GRB}}`$ for the total energy $`W_p`$ transferred to accelerated protons.

For the typical afterglow distance $`R_{\mathrm{sh}}10^{16}`$ cm, when the shock moves with constant Lorentz factor $`\mathrm{\Gamma }_{\mathrm{sh}}100`$, we shall use the following parameters for the upstream HE protons: the power-law spectrum $`E^{\gamma _g}`$ with the “standard” (see subsection VII.2.1) $`\gamma _g=2.2`$, with conservative $`E_{\mathrm{min}}\mathrm{\Gamma }_{\mathrm{sh}}^2E_i10^{13}`$ eV and with $`E_{\mathrm{max}}=eBR_{\mathrm{sh}}\mathrm{\Gamma }_{\mathrm{sh}}/\sqrt{3}1\times 10^{21}`$ eV, valid for very large $`B5`$ G. The fraction of the total energy $`W_p`$ in the form of particles with $`E_p1\times 10^{19}`$ eV is $`f_{\mathrm{UHECR}}=0.039`$. In exotic case of existence of pre-accelerated particles in pre-supernova wind with $`E_i10^{13}`$ eV, the minimal energy is $`E_{\mathrm{min}}10^{17}`$ eV and $`f_{\mathrm{UHECR}}=0.29`$.

As a matter of fact, part of accelerated particles are located downstream. Because they are isotropized in the rest frame, the minimum energy in the observer frame is shifted to

$$E_{\mathrm{min}}=\frac{m_pc^2}{2}\frac{\mathrm{\Gamma }_{\mathrm{sh}}}{\mathrm{\Gamma }_p^{}}(1+\mathrm{\Gamma }_p^{}_{}{}^{}2/\mathrm{\Gamma }_{\mathrm{sh}}^2),$$
(23)
where $`\mathrm{\Gamma }_p^{}`$ is the Lorentz factor of proton in the rest system. The energy (23) corresponds to a proton moving backward relative to shock. Additionally the number of particles with maximum energy diminishes due to angular distribution. It results in diminishing of $`f_{\mathrm{UHECR}}`$; we shall not take this effect into account.

There is one more effect which diminishes $`f_{\mathrm{UHECR}}`$ for downstream protons. The highest energy protons can scatter to backward direction and escape through the rear surface of the shell, while the low-energy protons are better confined in the downstream region. In the laboratory frame the backward moving protons have energy about 1 GeV, as given by Eq. (23). The difference of proton energy (in the laboratory frame) is transferred to the shell in the process of magnetic scattering. This energy will go back to particle acceleration, and most of it will be taken by low-energy particles.

Adiabatic energy losses, not taken into account in the relativistic shock acceleration, are discussed in previous subsection.

In realistic models the shock propagates in diminishing magnetic field, e.g. $`B(R)=B_0R_0/R`$ M. Vietri, D. De Marco and Guetta (2003), valid for the stellar wind. For such field it is easy to calculate the deflection angle $`\theta `$ for proton propagating from $`R_{\mathrm{sh}}`$ to $`R`$,

$$\theta \frac{eB_{\mathrm{sh}}R_{\mathrm{sh}}}{E}\mathrm{ln}R/R_{\mathrm{sh}},$$
(24)
and the deflection angle at which the shock overtakes the proton, $`\theta _c=\sqrt{3}/\mathrm{\Gamma }_{\mathrm{sh}}`$, is the same as for constant magnetic field. The maximum acceleration energy $`E_{\mathrm{max}}`$ determined from the condition $`\theta =\theta _c`$ remains also approximately the same. The exit of particles in the interstellar medium occurs when $`\mathrm{\Gamma }_{\mathrm{sh}}`$ diminishes, most notably at the Sedov phase.

In Fig. 20 the calculated spectrum of UHECR from acceleration by external shock is given by curve 2 for $`\gamma _g=2.2`$ and $`E_{\mathrm{max}}=1\times 10^{21}`$ eV. The UHECR emissivity needed to fit the flux at $`E=1\times 10^{19}`$ eV to the observational data is $`(>1\times 10^{19}\mathrm{eV})=4.7\times 10^{44}`$ erg Mpc-3 yr-1. The lower limit to HE proton emissivity (without adiabatic energy losses and other effects discussed above) is given by

$$_p>(>1\times 10^{19}\mathrm{eV})/f_{\mathrm{UHECR}}=1.2\times 10^{46}\mathrm{ergMpc}^3\mathrm{yr}^1$$
(25)
for $`E_{\mathrm{min}}=1\times 10^{13}`$ eV, and $`_p>1.6\times 10^{45}`$ erg Mpc-3 yr-1 for exotic case with $`E_{\mathrm{min}}=1\times 10^{17}`$ eV. The emissivity (25) is 2 - 3 orders of magnitudes higher than the observed GRB emissivity $`(0.613)\times 10^{43}`$ erg Mpc-3yr-1. The exotic assumption of $`E_{\mathrm{min}}=1\times 10^{17}`$ eV reduces this discrepancy only by one order of magnitude.

#### VII.3.3 Discussion and conclusion

The calculations above agree with earlier works V. Berezinsky, A. Z. Gazizov and Grigorieva (a); Stecker (2000); Scully and Stecker (2002); Stecker and Scully (2005) (note the emissivities given in the last work of this list).

The main difference between calculations V. Berezinsky, A. Z. Gazizov and Grigorieva (a); Stecker (2000); Scully and Stecker (2002); Stecker and Scully (2005)) on the one hand, and those by Waxman (e.g. Waxman (2001)) and Vietri (e.g. M. Vietri, D. De Marco and Guetta (2003)), on the other hand, consists in the following.

While Waxman and Vietri compare GRB emissivity with emissivity of the observed UHECR (Vietri at $`E1\times 10^{17}`$ eV), we compare it with emissivity, which corresponds to energy transferred to all protons in the process of their acceleration. Apart from spectral fraction of energy, it includes also the energy lost by protons before they escape to extragalactic space (adiabatic energy losses, energy transferred to the magnetized plasma outside the jet and others discussed above). Waxman for UHECR energy input uses the value $`E^2d\dot{n}_{\mathrm{CR}}/dE=0.8\times 10^{44}`$ erg Mpc-3yr-1 Waxman (2001), which after integration from $`E_{\mathrm{min}}=1\times 10^{11}`$ eV to $`E_{\mathrm{max}}=1\times 10^{21}`$ eV corresponds to emissivity $`_p=1.8\times 10^{45}`$ erg Mpc-3yr-1 to be compared with the observed GRB emissivity $`(0.613)\times 10^{43}`$ erg Mpc-3yr-1. Because of adiabatic energy losses, $`_p`$ is expected to be in fact much higher than the above-mentioned value.

Vietri *et al.* M. Vietri, D. De Marco and Guetta (2003) indicate UHECR energy release required by observations as $`1\times 10^{45}`$ erg Mpc-3yr-1 for energy interval $`1\times 10^{17}1\times 10^{21}`$ eV (compared with our value $`1.6\times 10^{45}`$ erg Mpc-3yr-1 for the same energy interval). If one includes the energy input to all accelerated protons, which accompany ultra-relativistic shock, from $`E_{min}\mathrm{\Gamma }_{\mathrm{sh}}^2E_i`$ to $`E_{\mathrm{max}}`$, this emissivity increases at least by factor 7.3, in conflict with observations by factor 60. There is some difference with our work in numerical calculations: Vietri *et al.* uses the HiRes data, we – energy-shifted data by AGASA and HiRes, which are in good intrinsic agreement (see Fig. 20); Vietri *et al.* introduce corrections due interaction fluctuations, we do not need to do it, because we normalize our calculations by observational data at $`E=2\times 10^{19}`$ eV, where fluctuations are absent. The maximum acceleration energy $`E_{\mathrm{max}}=1\times 10^{21}`$ eV at $`R_{\mathrm{sh}}=1\times 10^{16}`$ eV requires according to our calculations (it is by factor $`\sqrt{3}`$ less than in Gallant and Achterberg (1999)) magnetic field in the presupernova wind at this distance $`6`$ G, which is very large and marginally possible only for the SupraNova model Vietri and Stella (1998, 1999).

Our calculations require the UHECR emissivity $`(>1\times 10^{19}\mathrm{eV})=4.4\times 10^{44}`$ erg Mpc-3yr-1 for internal shock acceleration ($`\gamma _g=2.0`$ and $`E_{\mathrm{max}}=1\times 10^{21}`$ eV) and $`4.7\times 10^{44}`$ erg Mpc-3yr-1 for external ultra-relativistic shock acceleration ($`\gamma _g=2.2`$ and $`E_{\mathrm{max}}=1\times 10^{21}`$ eV). Including the energy injected in low energy protons (not necessarily coming outside), this emissivity should be a minimum of 10 - 30 times higher, and thus contradiction with GRB emissivity reaches conservatively the factor 30 - 100. A natural solution to this problem is given by an assumption that accelerated protons carry $`100`$ times more energy than electrons, whose radiation is responsible for GRBs. However, for $`4^{}`$ beaming it increases the required energy output of SN as GRB source from $`W_{\mathrm{SN}}5\times 10^{50}`$ erg to $`10^{52}10^{53}`$ erg, which creates the serious problem for SN models of GRBs. This problem exists even for the case when one takes into account only “observed” UHECR emissivity $`(>1\times 10^{19}\mathrm{eV})`$. Then the minimal energy balance for internal shock model includes $`_{\mathrm{GRB}}=1.3\times 10^{44}`$ erg Mpc-3yr-1, $`_{\mathrm{afterglow}}_{\mathrm{GRB}}`$, $`(>1\times 10^{19}\mathrm{eV})`$, and $`_\nu (>1\times 10^{19}\mathrm{eV})`$, which results in SN energy output $`4\times 10^{51}`$ erg. For external shock model the energy balance does not include neutrinos and it results in $`W_{\mathrm{SN}}3\times 10^{51}`$ erg.

## VIII  Conclusions

This work is naturally subdivided into two major parts: (i) the model-independent analysis of spectral signatures of UHE proton interaction with CMB, and (ii) the model-dependent analysis of transition from extragalactic to galactic cosmic rays and model-dependent study of most probable UHECR sources: AGN and GRBs.

In the first part the number of assumptions is minimal: we assume the power-law generation spectrum of extragalactic protons $`E^{\gamma _g}`$ and analyze uniform distribution of the sources with a possible distortion in the form of large-scale source inhomogeneities and local overdensity or deficit of the sources. The reference spectrum is given by the universal spectrum, calculated for the homogeneous distribution of the sources, when separation of sources $`d0`$. According to the propagation theorem Aloisio and Berezinsky (2004), in this case the spectrum does not depend on propagation mode, e.g. it is the same for rectilinear and diffusive propagation. We analyze also the spectra with different source separations, giving emphasis to separations $`d3050`$ Mpc, favorable by clustering of UHECR P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a); H. Yoshiguchi, S. Nagataki and Sato (2004a).

Calculation of the universal spectrum is based on equation of conservation of particles with continuous energy losses of protons interacting with CMB. We performed the new calculations of energy losses (first presented in V. Berezinsky, A. Z. Gazizov and Grigorieva (a)), which agree well with Berezinsky and Grigorieva (1988); T. Stanev *et al.* (2000), see Fig. 1. The technical element needed for calculations, the ratio of energy intervals at generation and observation $`dE_g/dE`$, is calculated by the new methods (see Appendix B). The analysis of spectral features is convenient to perform with help of modification factor, $`\eta (E)=J_p(E)/J_p^{\mathrm{unm}}(E)`$, where the spectrum $`J_p(E)`$ is calculated with all energy losses included, and unmodified spectrum $`J_p^{\mathrm{unm}}(E)`$ \- only with account of adiabatic energy losses (red shift).

There are four major signatures of UHE protons interacting with CMB: GZK cutoff, bump, dip and the second dip.

GZK cutoff Greisen (1966); Zatsepin and Kuzmin (1966) is the most prominent signature, which consists in the steepening of the proton diffuse spectrum due to photopion production. The beginning of the GZK cutoff at $`E_{\mathrm{GZK}}4\times 10^{19}`$ eV is difficult to observe. The quantitative characteristic of the GZK cutoff is given by energy $`E_{1/2}`$, where the flux with the cutoff becomes lower by factor 2 in comparison with power-law extrapolation from lower energies. This quantity is possible to observe only in integral spectrum. Practically independent of $`\gamma _g`$ $`E_{1/2}=5.3\times 10^{19}`$ eV. Data of Yakutsk and HiRes are compatible with this value (see Fig. 4), while data of AGASA contradict it. As to spectral shape of the GZK cutoff, we found it very model-dependent: it depends on separation of sources $`d`$, on local source overdensity or deficit, on maximum acceleration energy $`E_{\mathrm{max}}`$, and being so uncertain, the GZK steepening can be imitated by $`E_{\mathrm{max}}`$ acceleration steepening.

Bump is produced by pile-up protons, experiencing the GZK process. These protons do not disappear, they only shift in energy towards that where the GZK cutoff begins. The bump is very clearly seen in the spectra of individual sources (see left panel of Fig. 7), but they disappear in the diffuse spectrum (right panel of Fig. 7), because individual bumps are located at different energies.

Dip is the most robust signature of UHE protons interacting with CMB. It has very specific shape (see Fig. 2), which is difficult to imitate by other processes, unless one has a model with many free parameters. Dip is produced by pair production on the CMB photons $`p+\gamma p+e^++e^{}`$ and is located in energy interval $`1\times 10^{18}4\times 10^{19}`$ eV. Dip has two flattenings (see Fig. 2), one at energy $`E1\times 10^{19}`$ eV, which automatically reproduces the ankle, well seen in the observational data (see Fig. 8), and the other at $`E1\times 10^{18}`$ eV, which provides transition from extragalactic to galactic cosmic rays. The prediction of the dip shape is robust. It is not noticeably modified by many phenomena included in calculations: by different distances between sources (Fig. 9), by changing the rectilinear propagation of protons to the diffusive (Fig. 10), by different $`E_{\mathrm{max}}`$ (Fig. 5), by local overdensity and deficit of the sources (Fig. 6), by large-scale source inhomogeneities and fluctuations in $`p\gamma `$-interactions (Fig. 15).

The dip is very well confirmed by data of Akeno-AGASA, HiRes, Yakutsk and Fly’s Eye (see Fig. 19 for the latter), while data of Auger at this stage do not contradict the dip (Fig. 8). The agreement with the data of each experiment Akeno-AGASA, HiRes and Yakutsk is characterized by $`\chi ^21920`$ for about 20 energy bins, and with only two free parameters, $`\gamma _g`$ and the flux normalization constant. The best fit is reached at $`\gamma _g=2.7`$, with the allowed range 2.55 - 2.75.

The dip is used for energy calibration of the detectors, whose systematic energy errors are up to $`20\%`$. For each detector independently the energy is shifted by factor $`\lambda `$ to reach the minimum $`\chi ^2`$. We found $`\lambda _{\mathrm{Ag}}=0.9`$$`\lambda _{\mathrm{Hi}}=1.2`$ and $`\lambda _{\mathrm{Ya}}=0.75`$ for AGASA, HiRes and Yakutsk detectors, respectively. Remarkably, that after this energy shift the fluxes and spectra of all three detectors agree perfectly, with discrepancy between AGASA and HiRes at $`E>1\times 10^{20}`$ eV being not statistically significant (see Fig. 11). The AGASA excess over predicted GZK cutoff might have the statistical origin combined with systematic errors in energy determination, and the GZK cutoff may exist.

The difference in energy shifts between fluorescent energy measurement (HiRes) and on-ground measurements (AGASA, Yakutsk) with $`\lambda >1`$ for HiRes and $`\lambda <1`$ for AGASA and Yakutsk may signal a systematic difference in energy determination by these two methods.

The excellent agreement of the dip with observations should be considered as confirmation of UHE proton interaction with CMB.

For astrophysical sources of UHECR the presence of nuclei in primary radiation is unavoidable. Whatever a source is, the gas in it must contain at least helium of cosmological origin with ratio of densities $`n_{\mathrm{He}}/n_\mathrm{H}=0.079`$. Any acceleration mechanism operating in this gas must accelerate helium nuclei, if they are ionized. The fraction of helium in the primary radiation depends critically on the mechanism of injection. Fig. 13 shows that while 10% mixing of helium with hydrogen in the primary radiation leaves good agreement of dip with observations, 20% mixing upsets this agreement. One may conclude that good agreement of dip with observations indicates a particular acceleration mechanism, more specifically – an injection mechanism. One may also note that the shape of the dip contains the information about chemical composition of primary radiation (see Fig. 12).

The low-energy flattening of the dip is seen for both rectilinear and diffusive propagation (see Fig. 16). It is explained by transition from adiabatic energy losses to that due to $`e^+e^{}`$-pair production. Low-energy flattening of the extragalactic spectrum provides transition to more steep galactic component, as one can see it in the left panel of Fig. 17. Note, that when spectrum is multiplied to $`E^{2.5}`$ the flat spectrum looks like one raising with energy. In Fig. 17 the transition is shown for the diffusive spectrum. This model of transition implies that at energy below $`E_{\mathrm{tr}}5\times 10^{17}`$ eV the galactic iron is the dominant component, while at the higher energy extragalactic protons with some admixture of extragalactic helium dominate. The observational spectrum feature, which corresponds to this transition, is the second knee. The prediction of this model, the dominance of proton component at $`E>1\times 10^{18}`$ eV, is confirmed by data of HiRes Archbold and Sokolsky (2003); R. U. Abbasi *et al.* \[HiRes Collaboration\] (2005), HiRes-MIA T. Abu-Zaayad *et al.* \[HiRes Collaboration\] (2000), by Yakutsk A. V. Glushkov *et al.* \[Yakutsk collaboration\] (2000) and does not contradict to the Haverah Park data at $`E>(12)\times 10^{18}`$ eV M. Ave *et al.* \[Haverah Park collaboration\] (2003). However, the AkenoM. Honda *et al.* \[Akeno Collaboration\] (1993) and Fly’s Eye D. J. Bird *et al.* \[Fly’s Eye collaboration\] (1994) data favor the mixed composition dominated by heavy nuclei.

The second dip is produced by interplay between pair production and photopion production. It is a narrow feature at energy $`E_{2\mathrm{d}\mathrm{i}\mathrm{p}}=6.3\times 10^{19}`$ eV, which can be observed by detectors with good energy resolution. If detected, it gives the precise mark for energy calibration of a detector.

The UHECR sources have to satisfy two conditions: they must be very powerful and must accelerate particles to $`E_{\mathrm{max}}1\times 10^{21}`$ eV. There is one more restriction S. L. Dubovsky, P. G. Tinyakov and Tkachev (2000); Fodor and Katz (2001); H. Yoshiguchi, S. Nagataki and Sato (2003); G. Sigl, F. Miniati and Ensslin (2003a); H. Yoshiguchi, S. Nagataki and Sato (2004a); P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a), coming from observation of small-scale clustering: the space density of the sources should be $`(13)\times 10^5`$ Mpc-3, probably with noticeable uncertainty in this value. Thus, these sources are more rare, than typical representatives of AGN, the Seyfert galaxies, whose space density is $`3\times 10^4`$ Mpc-3. The sources could be the rare types of AGN, and indeed the analysis of Tinyakov and Tkachev (2001) show statistically significant correlation between directions of particles with energies $`(48)\times 10^{19}`$ eV and directions to AGN of the special type – BL Lacs. The acceleration in AGN can provide the maximum energy of acceleration to $`E_{\mathrm{max}}10^{21}`$ eV. The appropriate mechanisms can be acceleration by relativistic and ultra-relativistic shocks (see section VII.2.1). These mechanisms typically provide spectra flatter than one with $`\gamma _g2.62.8`$, needed for agreement of the dip with data. The distribution of the sources over $`E_{\mathrm{max}}`$ suggested in Kachelriess and Semikoz (2005a) solve this problem: starting from some energy $`E_c`$ the spectrum becomes steeper. In this case the diffuse spectrum below $`E_c`$ has canonical generation index $`\gamma _g=2.0`$ or $`\gamma _g=2.2`$, and above this energy $`\gamma _g2.62.8`$ to fit the observations. The energy $`E_c`$ is a free parameter of the models, and we keep it typically as $`10^{17}10^{18}`$ eV. Such broken generation spectrum relaxes the requirements for source luminosities. For example, for $`\gamma _g=2.7`$ and $`E_c=1\times 10^{18}`$ eV, the typical emissivity needed to fit the observed flux is $`_0=n_sL_p3\times 10^{46}`$ erg Mpc-3yr-1, which is held for the source density $`n_s=2\times 10^5`$ Mpc-3, estimated from the small-angle clustering P. Blasi and D. De Marco (2004); Kachelriess and Semikoz (2005a), and for the source luminosity $`L_p6\times 10^{43}`$ erg/s, very reasonable for AGN. The calculated spectra are in a good agreement with observations (see Fig. 19). The transition from galactic to extragalactic cosmic rays in this model occur at the second knee (see left panel of Fig. 17).

Another potential UHECR sources are GRBs. They have large energy output and can accelerate particles to UHEs. Acceleration occurs at the multiple inner shocks Waxman (1995, 2000) and at ultra-relativistic external shock Vietri (1995); M. Vietri, D. De Marco and Guetta (2003). The maximum acceleration energy can be as high as $`1\times 10^{21}`$ eV. In case of ultra-relativistic shock with Lorentz factor of shock $`\mathrm{\Gamma }_{\mathrm{sh}}100`$ and at distance $`R_{\mathrm{sh}}1\times 10^{16}`$ cm, typical for afterglow region, the maximum energy $`E_{\mathrm{max}}=(1/\sqrt{3})eBR_{\mathrm{sh}}\mathrm{\Gamma }_{\mathrm{sh}}`$ reaches $`1\times 10^{21}`$ eV in case of very strong external magnetic field $`B6`$ G. Such strong magnetic field can be marginally provided by two combined effects: strong field in presupernova wind (possible only in the model of SupraNova Vietri and Stella (1998, 1999)) combined with preshock amplification of the field.

The energy output of GRBs necessary to provide the observed UHECR flux imposes the serious problem Stecker (2000); Scully and Stecker (2002); Stecker and Scully (2005); V. Berezinsky, A. Z. Gazizov and Grigorieva (a). It is usually assumed that energy output in the form of accelerated protons is of order of magnitude of that observed in GRB photons. The observed emissivity of GRBs is

$$(0.613)\times 10^{43}\mathrm{ergMpc}^3\mathrm{yr}^1,$$
(26)
where the lower value is from works by Schmidt Schmidt (1999, 2001) and the large value - from work by Frail *et al.* D. A. Frail *et al.* (2001). Both authors used in their estimates the similar energy interval of photons, (10 - 1000 keV) and (0.2 - 2000 keV).

Our calculations of UHECR spectra for internal shocks (see Fig 20 curve 1) with $`\gamma _g=2.0`$ and $`E_{\mathrm{max}}=1\times 10^{21}`$ eV give UHECR emissivity $`(1\times 10^{19}\mathrm{eV})=4.4\times 10^{44}\mathrm{ergMpc}^3\mathrm{yr}^1`$, by factor 3.4 - 73 higher than (26). If to include $`_\nu _{\mathrm{UHECR}}`$ contradiction with (26) increases to 6.8 - 150. If to include in balance the protons with energies down to $`E_{\mathrm{min}}1\times 10^{11}`$ eV (not necessarily coming out) the contradiction increases by another factor 5. Adiabatic energy losses further increase this discrepancy.

Our calculations for external shock acceleration (see curve 2 in Fig 20) with $`\gamma _g=2.2`$ and $`E_{\mathrm{max}}=1\times 10^{21}`$ eV results in UHECR emissivity $`(1\times 10^{19}\mathrm{eV})=4.7\times 10^{44}\mathrm{ergMpc}^3\mathrm{yr}^1`$, i.e. by factor 3.6 - 78 times higher than (26). It increases by factor 26, if to include protons with energies down to $`E_{\mathrm{min}}\mathrm{\Gamma }_{\mathrm{sh}}^2E_i10^{13}`$ eV, or by factor 3.5 for unrealistic case of $`E_{\mathrm{min}}1\times 10^{17}`$ eV. Adiabatic energy losses are less than in case of internal shock acceleration.

This energy crisis can be resolved, assuming that accelerated protons in GRBs take away $`100`$ times more energy than electrons responsible for GRB photons. However, this assumption dramatically affects the status of SN origin for GRBs , even in the case of narrow $`4^{}`$ beaming.

### VIII.1 Predictions for the Auger measurements

What do we predict for the Auger observations on the basis of our analysis?

* We predict the dip at energy $`1\times 10^{18}4\times 10^{19}`$ eV not only as the spectrum shape, but with the absolute values of flux, based on the agreement of AGASA-HiRes-Yakutsk fluxes after energy calibration by the dip (see Fig. 11). This spectrum is displayed in the left panel of Fig. 21. The distortion of the shape of the predicted spectrum can occur only due to presence of extragalactic nuclei.
* The above-mentioned distortion is important as a tool for measuring the chemical composition of UHECR. The chemical composition is the most difficult problem for UHECR experiments, in which little progress has been reached during last 30 years. The measurement of the dip shape gives the precise method of determination of 10 - 15 % of admixture of helium in primary UHECR (see Fig. 12). The distortion of the proton dip most noticeably occurs at energy $`E>3\times 10^{19}`$ eV, (see Fig. 13) where precise measurements of Auger are expected.
* Normalization of energies measured by AGASA, Yakutsk and HiRes detectors with the help of the dip (see section IV.3) implies that energy measured by fluorescent method is higher than that by on-ground method (the energy shift for AGASA and Yakutsk is $`\lambda <1`$, while for HiRes $`\lambda >1`$). With help of the hybrid events Auger can resolve this problem.
* The second dip (see Fig. 24) may be marginally observed and used for energy calibration of the detector.
* Yakutsk and HiRes data give evidence for the numerical characteristic of the GZK cutoff $`E_{1/2}=5.3\times 10^{19}`$ eV in the integral spectrum (see Fig. 4). In the Auger data this value can be measured more reliably.
* The predicted shape of the GZK steepening has many uncertainties (see section III.1 for the discussion). However, beginning of the GZK cutoff in the narrow energy interval $`(47)\times 10^{19}`$ eV is predicted with the high accuracy (see left panel in Fig. 21). The precise spectrum measurement in this energy range can give the decisive proof of the GZK cutoff.
* The mass composition at $`1\times 10^{18}<E<1\times 10^{19}`$ eV measured by the fluorescent technique can distinguish between models of the second knee transition from galactic to extragalactic cosmic rays (protons with 10 - 20% of helium) and the ankle transition (iron dominance)
* The bumps in the spectrum of single sources, if detected, can be be observed (see Fig. 21).

###### Acknowledgements.

We acknowledge the participation of R. Aloisio and P. Blasi in some parts of this work. B. Hnatyk, D. De Marco and M. Kachelriess are thanked for useful discussions. Participation of M. I. Levchuk at the initial stage of pair-production energy loss recalculation is gratefully acknowledged. We thank ILIAS-TARI for access to the LNGS research infrastructure and for the financial support through EU contract RII3-CT-2004-506222. The work of S.G. is partly supported by grant LSS-5573.2006.2 .

## Appendix A  Calculations of energy losses

We consider here some details of calculations of energy-losses due to production of $`e^+e^{}`$-pairs and pions.

Pair production energy loss of ultrahigh-energy protons in low-temperature photon gas, e.g. CMB,

$$p+\gamma _{CMB}p+e^++e^{},$$
(27)
has been previously discussed in many papers. The differential cross-section for this process in the first Born approximation was originally calculated in 1934 by Bethe and Heitler Bethe and Heitler (1934) and Racah Racah (1934). In 1948 Feenberg and Primakoff Feenberg and Primakoff (1948) obtained the pair production energy loss rate using the extreme relativistic approximation for the differential cross-section. In 1970 the accurate calculation was performed by Blumenthal Blumenthal (1970). Later some analytical approximations to differential cross-sections were obtained in Ref. M. Chodorowski, A. A. Zdziarski and Sikora (1992).

All authors neglected the recoil energy of proton, putting $`m_p\mathrm{}`$, the effect being suppressed by a factor of $`m_e/m_p5\times 10^4`$.

In spite of the fact that all calculations actually used the same Blumenthal approach, there are noticeable discrepancies in the results of different authors; they are displayed in Fig. 1b of the Ref. T. Stanev, R. Engel, A. Muecke, R. J. Protheroe and Rachen (2000).

To clarify the situation we recalculate the pair production energy loss of high-energy proton in the low-temperature photon gas. In contrast to Ref. Blumenthal (1970) we use the first Born approximation approach of Ref. Berg and Linder (1961) taking into account the finite proton mass. The exact non-relativistic threshold formula with corrections due to different Coulomb interactions of electron and positron with the proton (see e.g. Ref. V. B. Berestetskii, E. M. Lifshits and Pitaevskii (1980)) is used. No series expansions of $`\sigma (\epsilon _r)`$, where $`\epsilon _r`$ is the photon energy in the proton rest system, are involved in our calculations.

The average fraction of proton energy loss in one collision with a photon is defined as

$$1x(\epsilon _r)=\frac{1}{\sigma (\epsilon _r)}\underset{x_{min}(\epsilon _r)}{\overset{x_{max}(\epsilon _r)}{}}𝑑x(1x)\frac{d\sigma (\epsilon _r,x)}{dx},$$
(28)
where $`x=E_p^{}/E_p`$, with $`E_p`$ and $`E_p^{}`$ being the incident and final proton energies, respectively, in the laboratory system. The upper and lower limits on fractions $`x`$ are given by

$$x_{_{\mathrm{min}}^{\mathrm{max}}}(\epsilon _r)=x_c(\epsilon _r)\pm y_c(\epsilon _r),$$
(29)
$$x_c(\epsilon _r)=\frac{1}{2}\left(1+\frac{m_p^2m_\pi ^2}{m_p^2+2m_p\epsilon _r}\right);y_c(E_c)=\sqrt{x_c^2(\epsilon _r)\frac{m_p^2}{m_p^2+2m_p\epsilon _r}}.$$
(30)
Our strategy was to calculate $`1x(\epsilon _r)`$ by performing the direct fourfold integration of the exact matrix element over the phase space. It should be noted, that direct numerical integration, especially at high energies, is difficult in this case because of forward-backward peaks in the electron-positron angular distributions. To overcome this problem, we performed two integrations over polar and azimuth angles in the $`e^+e^{}`$ subsystem analytically. This was facilitated by using of the *MATHEMATICA 4* code. The residual two integrations over energy and scattering angle in the initial $`p\gamma `$ subsystem were carried out numerically. We calculate simultaneously the total cross-section for pair production. The accuracy of our calculations was thus controlled by comparison of the calculated total cross-section with the well-known Bethe-Heitler cross-section.

The average fraction of proton energy lost in one collision with a photon is plotted in Fig. 22 as a function of the photon energy in the proton rest system.

The product of this fraction and the total cross-section for pair production is shown in Fig. 22. This function should be integrated over the photon spectrum to obtain the average energy loss due to pair production in the photon gas with this spectrum.

The comparison of our calculations with Ref. T. Stanev, R. Engel, A. Muecke, R. J. Protheroe and Rachen (2000) shows the negligible difference (see Fig. 1).

Calculations of photopion energy losses are described in section II.1. In the right panel of Fig. 22 we show $`1x\sigma `$ for these energy losses as function of $`\epsilon _r`$.

## Appendix B  Connection between energy intervals at epochs of production and observation

If we consider the protons with energy $`E`$ in the interval $`dE`$ at the epoch with red-shift $`z=0`$, what will be the corresponding interval of generation $`dE_g`$ at epoch $`z`$, when energy of a proton was $`E_g(z)`$? The connection between these two intervals is given by Eq. (36) of Ref. Berezinsky and Grigorieva (1988). Here we shall confirm this formula using a different, more simple derivation. Note, that intermediate formulae (40) and (41) used for derivation of final Eq. (36) in Ref. Berezinsky and Grigorieva (1988) have a misprint: the correct power of $`(1+z)`$ term there is 3, not 2.

Regarding the proton energy losses of a proton on CMB at arbitrary epoch with red-shift $`z`$, we shall use, as in Section II.1, the notation $`b(E,z)=dE/dt`$ and $`\beta (E,z)=(1/E)dE/dt`$ with $`b(E,z)=(1+z)^2b_0[(1+z)E]`$. Here and henceforth $`b_0(E)`$ and $`\beta _0(E)`$ are used for energy losses at $`z=0`$. The energy loss due to red-shift at the epoch $`z`$ and the Hubble parameter $`H(z)`$ are given by

$$\left(dE/dt\right)_{\mathrm{r}\mathrm{sh}}=EH(z),H(z)=H_0\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}.$$
(31)
The energy of a proton at epoch $`z`$,

$$E_g(z)=E+_t^{t_0}𝑑t\left[\left(dE/dt\right)_{\mathrm{r}\mathrm{sh}}+\left(dE/dt\right)_{\mathrm{CMB}}\right],$$
(32)
can be easily rearranged as

$$E_g(z)=E+_0^z\frac{dz^{}}{1+z^{}}E_g(z^{})+_0^z𝑑z^{}\frac{1+z^{}}{H(z^{})}b_0\left[(1+z^{})E_g(z^{})\right],$$
(33)
using $`dt/dz=H^1(z)/(1+z)`$ and Eq. (5).

Differentiating Eq. (33) with respect to $`E`$, one finds for energy interval dilation $`y(z)dE_g(z)/dE`$:

$$y(z)=1+_0^z\frac{dz^{}}{1+z^{}}y(z^{})+_0^z𝑑z^{}\frac{(1+z^{})^2}{H(z^{})}y(z^{})\left(\frac{db_0(E^{})}{dE^{}}\right)_{E^{}=(1+z^{})E_g(z^{})}.$$
(34)
Corresponding differential equation is

$$\frac{1}{y(z)}\frac{dy(z)}{dz}=\frac{1}{1+z}+\frac{(1+z)^2}{H(z)}\left(\frac{db_0(E^{})}{dE^{}}\right)_{E^{}=(1+z)E_g(z)}.$$
(35)
The solution of Eq. (35) is

$$y(z)\frac{dE_g(z)}{dE}=(1+z)\mathrm{exp}\left[\frac{1}{H_0}_0^z𝑑z^{}\frac{(1+z^{})^2}{\sqrt{\mathrm{\Omega }_m(1+z^{})^3+\mathrm{\Omega }_\mathrm{\Lambda }}}\left(\frac{db_0(E^{})}{dE^{}}\right)_{E^{}=(1+z^{})E_g(z^{})}\right],$$
(36)
where $`E_g(z)`$ is an energy at epoch z. In case of the highest energies the intergrand in Eq. (36) can be expanded in powers of $`z^{}`$, and after simple transformations one obtains

$$dE_g/dEb(E_g)/b(E),$$

valid with 10% accuracy at $`E>4\times 10^{19}`$ eV.

Eq. (36) coincides with Eq. (36) from Berezinsky and Grigorieva (1988).

We shall give now another derivation of Eq. (36) based on the exact solution to the kinetic equation for propagation of UHE protons with continuous energy losses. This solution automatically includes $`dE_g/dE`$ term which coincides with that given by Eq. (36). This equation for density of UHE protons, $`n_p(E,t)`$, reads

$$\frac{}{t}n_p(E,t)\frac{}{E}\left[\stackrel{~}{b}(E,t)n_p(E,t)\right]Q_g(E,t)=0,$$
(37)
with $`\stackrel{~}{b}(E,t)=EH(t)+b(E,t)`$, where $`EH(t)`$ describes adiabatic energy loss and $`b(E,t)`$ – those due interaction with CMB.

Written in the equivalent form

$$\frac{}{t}n_p(E,t)\stackrel{~}{b}(E,t)\frac{}{E}n_p(E,t)n_p(E,t)\frac{\stackrel{~}{b}(E,t)}{E}Q_g(E,t)=0,$$
(38)
this equation can be solved with help of auxiliary characteristic equation $`dE/dt=\stackrel{~}{b}(E,t)`$, which solution is $`(t)=E_g(E_i,t_i,t)`$, where like in Section II.2 $`E_g`$ is the generation energy at age $`t`$, if $`E=E_i`$ at $`t_i`$. When energy of a proton is taken on the characteristic $`E=(t)`$, the second term in the l.h.s. of Eq. (38) disappears and the solution found with help of integration factor is

$$n_p(,t)=_{t_{\mathrm{min}}}^t𝑑t^{}Q_g[(t^{}),t^{}]\mathrm{exp}\left[_t^{}^t𝑑t^{\prime \prime }\frac{}{}b((t^{\prime \prime }),t^{\prime \prime })\right].$$
(39)
Introducing the red-shift $`z`$ as a variable, $`dt/dz=H^1(z)/(1+z)`$, we obtain for present epoch ($`z=0`$):

$$n_p()=_0^{z_{\mathrm{max}}}\frac{dz^{}}{(1+z^{})H(z^{})}Q_g[(z^{}),z^{}]\mathrm{exp}\left[_0^z^{}\frac{dz^{\prime \prime }}{(1+z^{\prime \prime })H(z^{\prime \prime })}\frac{b(,z^{\prime \prime })}{}\right]$$
(40)
Using

$$\frac{b(,z)}{}=H(z)+(1+z)^3\left[\frac{b_0(E^{})}{E^{}}\right]_{E^{}=(1+z)(z)}$$

one finally gets

$$n_p()=\underset{0}{\overset{z_{\mathrm{max}}}{}}\frac{dz^{}}{(1+z^{})H(z^{})}Q_g[(z^{}),z^{}](1+z^{})\mathrm{exp}\left[\frac{1}{H_0}\underset{0}{\overset{z^{}}{}}𝑑z^{\prime \prime }\frac{(1+z^{\prime \prime })^2}{\sqrt{\mathrm{\Omega }_m(1+z^{\prime \prime })^3+\mathrm{\Omega }_\mathrm{\Lambda }}}\left(\frac{b_0(E^{})}{E^{}}\right)_{E^{}=(1+z^{\prime \prime })(z^{\prime \prime })}\right]$$
(41)
Since $`(z^{})=E_g(z^{})`$, one easily recognizes inside the integral the term $`dt^{}Q_g(E_g,t^{})`$ and $`dE_g/dE`$ as given by Eq. (36).

## Appendix C  Kinetic equation solution vs continuous loss approximation: physical interpretation

We shall compare here the kinetic equation solution for proton spectra, $`n_{kin}(E)`$, with that obtained in continuous-energy-loss (CEL) approximation, $`n_{cont}(E)`$, and study the physical meaning of their difference.

Additionally to the full kinetic equation (16) we shall consider here also the stationary kinetic equation with the photopion production only. Such equation describes realistically the spectrum at $`E>1\times 10^{20}`$ eV, where pion energy losses dominate, and effects of CMB evolution are negligible. However, we formally consider this equation at the lower energies, too, to study the difference between kinetic equation solution and CEL approximation.

The general kinetic equation is reduced to the stationary one with photopion interaction only, putting in Eq. (16) $`n/t=0`$, $`H(t)=0`$, $`T=T_0`$ and $`b_{pair}=0`$.

Then the stationary kinetic equation reads

$$P(E)n_p(E)+\underset{E}{\overset{E_{max}}{}}𝑑E^{}P(E^{},E)n_p(E^{})+Q_{gen}(E)=0,$$
(42)
with $`P(E)`$ and $`P(E^{},E)`$ given by Eqs. (17) and (18). Introducing variable $`x=E/E^{}`$, we rearrange the regeneration term into

$$\frac{1}{E}\underset{0}{\overset{1}{}}𝑑x\left(\frac{E}{x}\right)P(E/x,x)n_p(E/x).$$
(43)
Expanding the integrand into the Taylor series in powers of $`(1x)`$ one observes the exact compensation of the zero power term with $`P(E)n(E)`$. The first power term results in the stationary continuous energy loss equation

$$\frac{}{E}\left[b(E)n_p(E)\right]+Q_{gen}(E)=0,$$
(44)
and account of $`(1x)^2`$ term gives the stationary Fokker-Planck equation

$$\frac{}{E}\left[b(E)n_p(E)+\frac{}{E}\left(E^2D(E)n_p(E)\right)\right]+Q_{gen}(E)=0,$$
(45)
where the term with diffusion coefficient $`D(E)`$ in momentum space describes fluctuations.

This result reproduces the general proof (see e.g. Lifshitz and Pitaevskii (1981), p. 89) for the case when all processes go with a small energy transfer.

Let us now consider the analytic solution to Eq. (42) with $`E_{\mathrm{max}}\mathrm{}`$ in the asymptotic limit $`EE_{GZK}`$. We shall use here (and only here) a toy model with simplified assumptions about interactions, namely with one-pion production $`p+\gamma \pi +N`$, with isotropic pion distribution in c.m. system and with asymptotic cross-section $`\sigma _{p\gamma }\text{const}`$, when $`E_c\mathrm{}`$. With these assumptions we readily obtain

$$P(E)c\sigma _{p\gamma }n_\gamma ,P(E^{},E)c\sigma _{p\gamma }n_\gamma /E^{},$$
(46)
and using $`Q_{gen}(E)=Q_0E^{\gamma _g},`$ we arrive at

$$n_{kin}(E)=\frac{\gamma _g}{\gamma _g1}\frac{1}{\sigma _{p\gamma }n_\gamma c}Q_0E^{\gamma _g}.$$
(47)
For continuous energy loss equation (44), using $`b(E^{})=E^{}Ec\sigma _{p\gamma }n_\gamma `$ and $`E^{}E=\frac{1}{2}E^{}`$, valid for isotropic pion distribution in c.m. system, we obtain

$$n_{cont}(E)=\frac{2}{c\sigma _{p\gamma }n_\gamma }Q_0E^{\gamma _g}.$$
(48)
Therefore, the asymptotic ratio is

$$\frac{n_{kin}(E)}{n_{cont}(E)}=\frac{\gamma _g}{2},$$
(49)
valid also for $`\gamma _g=2`$. This result has been obtained earlier in V. S. Berezinsky, S. I. Grigorieva and Zatsepin (1975). Of course, at the considered energies nothing but fluctuations are responsible for ratio (49). This ratio, being asymptotic, can be used, however, as an estimate of role of fluctuations at $`E>>E_{GZK}`$.

Let us now come over to the case of finite $`E_{max}`$ and consider $`E`$ very close to $`E_{max}`$.

The regeneration term in Eq. (42) tends to zero and

$$n_{kin}(E)\frac{Q_0E_{max}^{\gamma _g}}{P(E_{max})},$$
(50)
having non-zero value at $`E=E_{max}ϵ`$ and being zero at $`E=E_{max}+ϵ`$, where $`ϵ`$ is infinitesimally small energy.

For particle density at $`E=E_{max}ϵ`$ in continuous loss approximation one obtains from Eq. (44)

$$n_{cont}(E_{max})=\frac{1}{b(E_{max})}_{E_{max}E}^{E_{max}}Q_0E^{\gamma _g}𝑑E\frac{ϵ}{E_{max}}0.$$
(51)
Therefore, $`n_{kin}(E)/n_{cont}(E)\mathrm{}`$ as $`EE_{max}`$.

This result has a clear physical meaning. In our model sources are homogeneously distributed and from nearby sources the particles with energy $`E_{max}`$ can arrive at the observational point with the same energy due to fluctuations. But in CEL approximation they cannot do it, since they loose energy at any distance traversed, whatever small it is. The described effect influences the ratio $`r(E)=n_{kin}(E)/n_{cont}(E)`$ at all energies close enough to $`E_{max}`$.

Let us now come over to the solution of Eq. (42) with interaction given by the detailed calculations as described in Section V. The calculated ratio $`r(E)=n_{kin}(E)/n_{cont}(E)`$ is plotted in Fig. 23 for $`E_{\mathrm{max}}=1\times 10^{23}`$ eV and for three values of $`\gamma _g`$ equal to 2.3, 2.7 and 2.9. One may observe in Fig. 23 some features described by our toy model: the ratio $`r(E)`$ increases with $`\gamma _g`$ as it should according to Eq. (49) and it tends to $`\mathrm{}`$ at $`E_{max}`$. For spectra of interest for the present work with $`2.3<\gamma _g<2.8`$ the ratio $`r(E)`$ at $`1\times 10^{19}E1\times 10^{21}`$ eV does not exceed $`10`$ %. Note that in Fig. 15 it is smaller at low energies because the non-fluctuating pair-production process is included there.

The ratio of the solutions for kinetic equation and for CEL equation includes fluctuations, but it includes also some other effects. The kinetic equation tends to CEL equation (44) in the limit of small $`1x`$. But in fact for all reasonable energies $`1x`$ is limited by minimum value $`m_\pi /(m_\pi +m_p)0.13`$, which cannot be considered as small value for the discussed accuracy. The same argument is valid for comparison of Monte-Carlo simulation with CEL approximation: why they should coincide if the energy transfer process is discrete? However, we do see from Fig. 23 that the ratio $`r(E)`$ differs from 1 just by a few percent, when energy becomes less than $`1\times 10^{19}`$ eV.

Appendix D gives an example when kin/CEL ratio at energy close to $`1\times 10^{20}`$ eV is explained not by fluctuations. At energy $`E6\times 10^{19}`$ eV the kin/CEL ratio is affected by an interplay between pair- and photopion production.

We shall now discuss the comparison of our calculations with Monte-Carlo (MC) simulations.

As was emphasized in Section V the kinetic equation solution must coincide exactly with MC simulation when the number of the MC tests tends to infinity, the distribution of sources is homogeneous, $`\gamma _g`$ and $`E_{\mathrm{max}}`$ are the same and interactions are included in the identical way. D. De Marco and M. Kachelriess run their MC programs specially for comparison with our calculations using $`\gamma _g=2.7`$, $`E_{\mathrm{max}}=1\times 10^{23}`$ eV and homogeneous distribution of sources. They used SOPHIA interaction model T. Stanev, R. Engel, A. Muecke, R. J. Protheroe and Rachen (2000), which is not identical to our interaction model especially at the highest energies $`E1\times 10^{21}`$ eV, but the average energy losses coincide well in both models (see Fig. 1). The comparison of MC/CEL and kin/CEL ratios are given in the Table below. The data of MC are kindly provided by D. De Marco (D) and M. Kachelriess (K).

From Table 2 one can see the tendency of smaller kin/CEL ratios in comparison with MC/CEL, which can be due to different interaction models used in these calculations. We are planning the joint work with D. De Marco and M. Kachelriess on the detailed study of the discussed effects by methods of kinetic equation and Monte-Carlo.

The main conclusion of our study is that fluctuations in photopion production modify only weakly the spectra with different $`\gamma _g`$ at energies up to $`1\times 10^{21}`$ eV (see Fig. 23).

## Appendix D  The second dip

As explained qualitatively in Section III.4 the second dip appears due to breaking of compensation between particle exit and regeneration given by the first and second terms in Eq. (42). This is a narrow feature near the energy $`E_{\mathrm{eq2}}`$ which is produced because of the sharp (exponential) increase of photopion energy losses with energy, seen in Fig. 1. At energy $`E`$ just slightly below $`E_{\mathrm{eq2}}`$ the continuous (adiabatic and pair-production) energy losses dominate, and spectrum is universal. At energy slightly above $`E_{\mathrm{eq2}}`$ the continuous energy losses are small, and spectrum is determined by Eq. (42) with the high accuracy compensation between the absorption, $`P(E)n_p(E)`$, and regeneration terms. At $`EE_{\mathrm{eq2}}`$ the continuous energy losses break this compensation, increasing the absorption term, and triggering thus the appearance of the dip. We shall study here this triggering mechanism quantitatively.

Coming back to the basic kinetic equation (16), we rearrange the terms connected with expansion of the universe (proportional to $`H(t)`$) and with pair-production energy losses, including them in $`P_{eff}(E,t)`$. We obtain

$$\frac{n_p(E,t)}{t}=P_{eff}(E,t)n_p(E,t)+\underset{E}{\overset{E_{max}}{}}𝑑E^{}P(E^{},E,t)n_p(E^{},t)+Q_{gen}(E),$$
(52)
where $`P_{eff}(E,t)=P(E,t)+P_{cont}(E,t)`$, with $`P(E,t)`$ given as before by Eq. (17) and

$$P_{cont}(E,t)=2H(t)\left[\beta _{pair}(E,t)+H(t)\right]\frac{\mathrm{ln}n(E,t)}{\mathrm{ln}E}\frac{b(E,t)}{E}$$
(53)
As was described above, $`P_{cont}(E)`$ breaks the compensation between $`P(E)n_p(E)`$ and the regeneration term in Eq. (52), triggering thus the modification of $`n_{kin}(E)`$. It is convenient to introduce the auxiliary trigger function defined at $`t=t_0`$ as

$$T(E)=\{\begin{array}{cc}P_{eff}(E)/P_{cont}(E)\hfill & \mathrm{at}EE_c\hfill \\ P_{eff}(E)/P(E)\hfill & \mathrm{at}EE_c,\hfill \end{array}$$
(54)
with $`E_c=6.1\times 10^{19}`$ eV, determined from equation $`P(E)=P_{cont}(E)`$ and being equal to $`E_{\mathrm{eq2}}`$.

The trigger function describes how $`P_{eff}(E)`$ in Eq. (52) is changing from $`P_{cont}(E)`$ at $`EE_c`$, where $`T=1`$, to $`P(E)`$ at $`EE_c`$, where $`T=1`$ as well. The calculated trigger function is plotted in the left panel of Fig. 24. In the calculations we used for $`n_p(E)`$ the universal spectrum, because as Fig. 15 shows, its distortion is small. The trigger function $`T`$ reaches 1 at $`E<3\times 10^{19}`$ eV, and therefore at these energies $`P_{eff}(E)P_{cont}(E)`$; the regeneration term is small at these energies, too, and thus $`n_{kin}(E)n_{cont}(E)`$ holds. At $`E>1\times 10^{20}`$ eV $`T1`$, $`P_{eff}(E)P(E)`$ and $`n_{kin}(E)`$ becomes the solution of Eq. (42) with photopion energy losses only and with $`n_{kin}(E)/n_{cont}(E)`$ as shown in Fig. 23. The maximum deflection of $`n(E)`$ from these two regimes reached at $`E_c=6.1\times 10^{19}`$ eV. For the steep spectrum $`n(E)`$, which is the realistic case, the second term in Eq. (53) is positive and large in comparison with other terms, $`P_{cont}(E)>0`$, the effective exit probability $`P(E)+P_{cont}(E)`$ increases, and thus $`n_{kin}(E)`$ decreases. The triggering mechanism predicts that $`E_c`$ does not depend on $`\gamma _g`$ and $`E_{\mathrm{max}}`$.

The exact calculations for $`n_{kin}(E)/n_{cont}(E)`$ are shown in the right panel of Fig. 24. In fact, this is just the ratio of the modification factors from Fig. 15 (’fluct.’ and ’cont.loss.’ curves), shown with a strong multiplication in order to demonstrate the fine structure of this function. One can observe the narrow dip at $`E_{2\mathrm{d}\mathrm{i}\mathrm{p}}=6.3\times 10^{19}`$ eV, which is only marginally seen in Fig. 15. The second dip appears in the solution of the Fokker-Planck equation, too.

The minimum of the dip coincides well with the maximum of the trigger function $`E_c=6.1\times 10^{19}`$ eV, and the small difference between them is caused by the fact that $`E_c`$ is calculated for $`t=t_0`$, while for $`E_{2\mathrm{d}\mathrm{i}\mathrm{p}}`$ the contribution from epochs with $`t<t_0`$ is present. The widths of both features are numerically similar. As the exact calculations show, the position of the dip does not depend on $`\gamma _g`$ and $`E_{\mathrm{max}}`$, in accordance with triggering mechanism.

We would like to emphasize the similarity of the first and the second dip. The first dip starts at energy $`E_{\mathrm{eq1}}`$, where pair-production energy losses become equal to the adiabatic energy losses. The second dip occurs near energy $`E_{\mathrm{eq2}}`$, where pair-production and pion-production energy losses are equal. Both dips are the faint features, not seen well when measured spectrum is exposed in a natural way $`J_{\mathrm{obs}}(E)`$ vs $`E`$. The first dip is clearly seen in ratio of measured spectrum to unmodified spectrum, $`J_{\mathrm{obs}}(E)/J_{\mathrm{unm}}(E)`$, the second dip is predicted to be visible in the ratio of the measured spectrum to the smooth universal spectrum, $`J_{\mathrm{obs}}(E)/J_{\mathrm{univ}}(E)`$. The latter is predicted as ratio $`n_{kin}(E)/n_{cont}(E)`$ shown in Fig. 24.

## Appendix E  The generation function for case of distribution of sources over maximal energies

We calculate here analytically the generation function $`Q_{gen}(E)`$, which gives the number of produced particles per unit comoving volume of universe and per unit time at $`z=0`$, in case the sources are distributed over maximum of acceleration $`\epsilon E_{max}`$.

We consider first the acceleration by non-relativistic shock, assuming the power-law spectrum of accelerated particles $`E^2`$ with exponential suppression $`\mathrm{exp}(E/\epsilon )`$ at $`E\epsilon `$. Namely, we take the generation function for a single source, $`q_{gen}(E)`$, which gives the number of particles with energy $`E`$ produced per unit time, in the form

$$q_{gen}(E,\epsilon )=\frac{L_p(\epsilon )}{2+\mathrm{ln}(\epsilon /E_{min})}\phi _{gen}(E,\epsilon ),$$
(55)
where $`E_{min}1`$ GeV is the minimum energy in acceleration spectrum, $`L_p(\epsilon )`$ is the luminosity of a source in the form of accelerated particles with energy $`E`$, and

$$\phi _{gen}(E,\epsilon )=\{\begin{array}{cc}E^2\hfill & \text{ at }E\epsilon ,\hfill \\ \epsilon ^2\mathrm{exp}(1E/\epsilon )\hfill & \text{ at }E\epsilon .\hfill \end{array}$$
(56)
Then the generation function per unit volume $`Q_g(E)`$ is

$$Q_{gen}(E)=_{\epsilon _{min}}^{\epsilon _{max}}𝑑\epsilon \frac{n_s(\epsilon )L_p(\epsilon )}{2+\mathrm{ln}(\epsilon /E_{min})}\phi _{gen}(E,\epsilon ),$$
(57)
where $`n_s(\epsilon )n_s[\epsilon ,L_p(\epsilon )]`$ is the space density of the sources.

We introduce also the spectral emissivity, $`(\epsilon )=n_s(\epsilon )L_p(\epsilon )`$ and the total emissivity $`_0=(\epsilon )𝑑\epsilon `$.

$`Q_{gen}(E)`$ can be evaluated from Eqs. (57) and (56) as follows:
for $`E<\epsilon _{min}`$,

$$Q_{gen}(E)E^2;$$
(58)
for $`\epsilon _{min}E\epsilon _{max}`$ the main contribution to $`Q_{gen}(E)`$ is given by

$$Q_{gen}^{(1)}(E)E^2_E^{\epsilon _{max}}𝑑\epsilon (\epsilon ),$$
(59)
while the integral from $`\epsilon _{min}`$ to $`E`$ results in

$$Q_{gen}^{(2)}(E)E^1(E);$$
(60)
and for $`E>\epsilon _{max}`$

$$Q_{gen}(E)\mathrm{exp}(E/\epsilon _{max}).$$
(61)
The physically interesting regime is given by $`\epsilon _{min}<E<\epsilon _{max}`$. To reproduce $`Q_{gen}(E)E^{2.7}`$ needed for the best fit of the dip, $`(\epsilon )`$ from (59) and (60) must be power-law function $`E^\beta `$ with $`\beta =1.7`$.

In this case we obtain

$$Q_{gen}(E)\{\begin{array}{cc}E^2\hfill & \text{ at }E<\epsilon _{min}\hfill \\ E^{2.7}\hfill & \text{ at }\epsilon _{min}E\epsilon _{max}\hfill \\ \mathrm{exp}(E/\epsilon _{max})\hfill & \text{ at }E>\epsilon _{max},\hfill \end{array}$$
(62)
which coincides with the phenomenological broken generation spectrum (20).

The function $`Q_{gen}(E)`$ from (62) can be easily normalized using the condition $`EQ_{gen}(E)𝑑E=_0`$.

One may generalize this derivation for the case

$$\phi _{gen}(E,\epsilon )=\{\begin{array}{cc}E^{(2+\alpha )}\hfill & \text{ at }E<\epsilon \hfill \\ E^{(2+\alpha )}f(\frac{E}{\epsilon })\hfill & \text{ at }E\epsilon ,\hfill \end{array}$$
(63)
where $`f(x)`$ is a function which diminishes with $`x`$ like $`e^x`$ or faster and $`f(1)=1`$.

Similar to (55) and (57) we have

$`q_{gen}(E,\epsilon )`$ $`=`$ $`\alpha E_{min}^\alpha L_p(\epsilon )\phi _{gen}(E,\epsilon ),`$ (64)
$`Q_{gen}(E,\epsilon )`$ $`=`$ $`\alpha E_{min}^\alpha {\displaystyle _{\epsilon _{min}}^{\epsilon _{max}}}𝑑\epsilon (\epsilon )\phi _{gen}(E,\epsilon ).`$ (65)
Using the properties of $`f(x)`$ we can prove again that the power-law function $`Q_{gen}(E)E^{\gamma _g}`$ in the interval $`\epsilon _{min}<E<\epsilon _{max}`$ needs the power-law function $`(\epsilon )\epsilon ^\beta `$.

Then using again the properties of $`f(x)`$ we obtain

$$Q_{gen}(E)\{\begin{array}{cc}E^{(2+\alpha )}\hfill & \text{ at }E\epsilon _{min}\hfill \\ E^{(1+\alpha +\beta )}\hfill & \text{ at }\epsilon _{min}E\epsilon _{max}\hfill \\ f(E/\epsilon _{max})\hfill & \text{ at }E>\epsilon _{max},\hfill \end{array}$$
(66)
We are interested in particular in case of acceleration at relativistic shock, when $`\alpha 0.20.3`$. The case of $`\gamma _g=1+\alpha +\beta =2.7`$ results then in $`\beta =1.7\alpha 1.41.5`$.