Acta Polytechnica doi:10.14311/AP.2013.53.0703 Acta Polytechnica 53(Supplement):703–706, 2013 © Czech Technical University in Prague, 2013 available online at http://ojs.cvut.cz/ojs/index.php/ap COMPUTATIONAL SCHEMES FOR THE PROPAGATION OF ULTRA HIGH ENERGY COSMIC RAYS Roberto Aloisioa,b,∗ a INAF, Osservatorio Astrofisico Arcetri, I-50125 Arcetri (Firenze) Italy b INFN, Laboratori Nazionali del Gran Sasso, I-67010 Assergi (L’Aquila) Italy ∗ corresponding author: aloisio@arcetri.astro.it Abstract. We discuss the problem of ultra high energy particles propagation in astrophysical backgrounds. We present two different computational schemes based on kinetic and Monte Carlo approaches. The kinetic approach is an analytical computation scheme based on the hypothesis of continuos energy losses while the Monte Carlo scheme takes into account also the stochastic nature of particle interactions. These schemes, which give quite reliable results, enable the computation of fluxes keeping track of the different primary and secondary components, providing a fast and useful workbench for studying Ultra High Energy Cosmic Rays. Keywords: particles astrophysics, ultra high energy cosmic rays, astrophysical backgrounds. 1. Introduction Ultra High Energy Cosmic Rays (UHECR) are the most energetic particles observed in nature, with en- ergies up to a several 1020 eV. Experimental studies of UHECR are currently being conducted in three different experiments: Auger in Argentina, HiRes and Telescope Array in the USA. The propagation of UHECR from the source to the observer is conditioned by their interactions with astro- physical backgrounds: the Cosmic Microwave Back- ground (CMB) and the Extragalactic Background Light (EBL). Understanding the key features of prop- agation is of paramount importance for interpreting experimental observations paving the way for the dis- covery of the astrophysical origin of these fascinating particles. Several features of the observed spectrum can be linked directly to the chemical composition of UHECR and to their sources [1, 2, 7, 9, 13]. One of the par- ticularly important features is the Greisin, Zatsepin and Kuzmin (GZK) suppression of the flux, an abrupt depletion of the observed proton spectrum, arising at energies E ' 5 × 1019 eV, due to the interaction of UHE protons with the CMB radiation field [9, 13]. GZK suppression, as follows from the original pa- pers, refers to protons and it is due to the photo- pion production process on the CMB radiation field (p+γCMB → π+p). In the case of nuclei the expected flux also shows a suppression at the highest energies which, depending on the nucleus species, is due to the photo-disintegration process on the CMB and EBL radiation fields (A+γCMB,EBL → (A−nN) +nN) [3]. Another important feature in the spectrum that can be directly linked with the nature of the primary par- ticles and their origin (galactic/extra-galactic) is the pair-production dip [1, 7]. This feature is present only in the spectrum of UHE extragalactic protons and, like GZK, is a direct consequence of the interaction of the proton with the CMB radiation field. In particular the dip brings a direct imprint of the pair production process p + γCMB → p + e + + e− suffered by protons. From the experimental point of view the situation is far from being clear with different experiments claiming contradictory results. The HiRes experi- ment, which is no longer taking data, showed a proton dominated spectrum till the highest energies [10, 11] while the Auger observations show a heavy mass com- position at energies E > 4 × 1018 eV [5]. This puzzling situation, with different experiments favoring different scenarios, shows once more the im- portance of a systematic study of UHECR propagation in astrophysical backgrounds. In the present paper we will review the main points of two alternative compu- tation schemes which enable the determination of the fluxes expected on earth fixing the injection spectrum and the distribution of sources. These two schemes are based on different approaches to modeling the interactions between particles and backgrounds: the continuum energy losses (CEL) approximation, which forms the basis of the kinetic approach, and the Monte Carlo (MC) technique. As we will discuss in the following these two different schemes give reliable results that, in the framework of different assumptions, agree each other and offer a suitable theoretical framework to study experimental results unveiling the intimate nature of UHECR. 2. Kinetic Equations The main assumption under which the kinetic the- ory is built is the CEL approximation [6], through which particle interactions are treated as a continuum process that continuously depletes the energy of the particles. 703 http://dx.doi.org/10.14311/AP.2013.53.0703 http://ojs.cvut.cz/ojs/index.php/ap Roberto Aloisio Acta Polytechnica UHECR propagating through astrophysical back- grounds suffer different interaction processes: • protons – UHE protons interact only with the CMB radiation field giving rise to the two processes of pair production and photo-pion production. Both of these reactions can be treated in the CEL hy- pothesis. • nuclei – UHE nuclei interact with the CMB and EBL radiation fields, suffering the process of pair production, for which only CMB is relevant, and photo-disintegration, which involves both CMB and EBL backgrounds. While the first process can be treated in the CEL hypothesis, the nucleus species being conserved, the second cannot be, producing a change in the nucleus species. Following Aloisio et al. [3], in the framework of the kinetic approach, we will treat the photo-disintegration process as a “decaying” process that simply depletes the flux of the propagating nucleus. Taking into account all energy loss processes we can describe the propagation of protons and nuclei through kinetic equations of the type: ∂np(Γ,t) ∂t − ∂ ∂Γ [bp(Γ,t)np(Γ,t)] = Qp(Γ,t) (1) ∂nA(Γ,t) ∂t − ∂ ∂Γ [nA(Γ,t)bA(Γ,t)]+ nA(Γ,t) τA(Γ,t) = QA(Γ,t) (2) where n is the equilibrium distribution of particles, b are the energy losses (adiabatic expansion of the Universe and pair/photo-pion production for protons or only pair-production for nuclei) Q is the injection of freshly accelerated particles and, in the case of nuclei, also the injection of secondary particles produced by photo-disintegration (see below). The energy losses b for protons or nuclei depend only on the CMB field and in the CEL hypothesis they can be computed analytically [1, 3, 7]. The second process that affects nuclei propagation is photo-disintegration over CMB and EBL backgrounds. This process is treated as a decaying process that depletes the flux of nuclei. It enters in the kinetic equation (see Eq. 2) through a sort of “life-time” of the nucleus under the photo-disintegration process. This “life-time” corresponds to the mean time needed for a nucleus of Lorentz factor Γ and atomic mass number A to lose, at least, one of its nucleons: 1 τA = c 2Γ 2 ∫ ∞ �0(A) d�rσ(�r,A)ν(�r)�r ∫ ∞ �r/(2Γ) d� nbkg(�) �2 (3) where σ(�r,A) is the photo-disintegration cross-section and ν(�r) is the multiplicity associated with this process, namely the average number of nucleons ex- tracted from the nucleus by a single interaction and nbkg = nCMB + nEBL. The dependence on red-shift of τA follows directly from the evolution with red-shift of the background photon densities nCMB and nEBL. In the case of CMB this dependence is known ana- lytically while for EBL one should refer to evolution models (in our computations we have used the model by Stecker et al. [12]). One important feature of the photo-disintegration process is that it starts to contribute to the propa- gation of nuclei at a Lorentz factor that is almost independent of the nuclei species Γcr ' 2 × 109 [3]. This is an important general characteristic of nuclei photo-disintegration process from which we can imme- diately deduce the dependence on the nuclei species of the energy corresponding to the photo-disintegration suppression of the flux: EAcut = AmNΓcr. A being the atomic mass number of the nucleus and mN the proton mass. From this expression for EAcut it is evident how the flux behavior could provide informations on the chemical composition of UHECR. In the case of He- lium (A = 4), suppression is expected around energies E ' 1019 eV while in the case of Iron (A = 56) sup- pression is expected at higher energies E ' 1020 eV. Let us now discuss the generation function QA(Γ,t) on the right hand side of Eq. 2. One should distin- guish between primary nuclei, i.e. nuclei accelerated at the source and injected in the intergalactic space, and secondary nuclei and nucleons, i.e. particles pro- duced as secondaries in the photo-disintegration chain. In the case of primaries the injection function is an assumption of the source model, while the injection of secondaries should be modeled taking into account the characteristics of the photo-disintegration pro- cess. The dominant process of photo-disintegration is one nucleon (N) emission, namely the process (A + 1) + γbkg → A + N. This follows directly from the behavior of the photo-disintegration cross-section (see [3, and references therein]) which shows the giant dipole resonance corresponding to one nucleon emis- sion. Moreover, at the typical energies of UHECR (E > 1017 eV) one can safely neglect the nucleus re- coil so that photo-disintegration will conserve the Lorentz factor of the particles. The production rate of secondary A-nucleus and A-associated nucleons will therefore be given by QA(Γ,z) = QAp (Γ,z) = nA+1(Γ,z) τA+1(Γ,z) (4) where τA+1 is the photo-disintegration life-time of the nucleus father (A + 1) and nA+1 is its equilib- rium distribution, the solution of the kinetic equation (Eq. 2). Using Eq. 4 we can build a system of coupled differential equations that starting from the pri- mary injected nuclei (A0) follows the complete photo- disintegration chain for all secondary nuclei (A < A0) and nucleons. Clearly secondary proton1 propagation will be described by the proper kinetic equation (Eq. 1) with an injection term given by Eq. 4. The solution 1Neutrons decay very fast into protons, so we will always refer to secondary protons. 704 vol. 53 supplement/2013 Computational Schemes for the Propagation of Ultra High Energy Cosmic Rays (E/eV) 10 log 18 18.5 19 19.5 20 20.5 ] -1 y r -1 sr -2 k m 2 J [e V 3 E 33 10 34 10 35 10 36 10 37 10 56 50 40 30 20 10 Figure 1. Flux of iron and secondary nuclei (A = 50, 40, 30, 20, 10) at z = 0 in the case of pure iron injection at the source with a power law injection index γ = 2.2. Full squares correspond to the SimProp result [4] while continuous lines correspond to the solution of the nuclei kinetic equation of [3]. of the kinetic equation for protons and nuclei can be worked out analytically. In the case of protons: np(Γ,z) = ∫ zmax z dz′ (1 + z′)H(z′) Qp(Γ ′,z) dΓ ′ dΓ , (5) being Qp the injection of primary protons or secondary protons (Eq. 4) and Γ ′ = Γ ′(Γ,z) is the characteristic function of the kinetic equation [3]. In the case of nuclei: nA(Γ,z) = = ∫ zmax z dz′ (1 + z′)H(z′) QA(Γ ′,z) dΓ ′ dΓ e−ηA(Γ ′,z′), (6) being, again, QA the injection of primary or secondary (Eq. 4) nuclei. The exponential term in Eq. 6 repre- sents the survival probability during the propagation time t′ − t for a nucleus with fixed A and can be computed according to Aloisio et al. [3]. The deriva- tive term dΓ ′/dΓ present in both solutions Eq. 5 and Eq. 6 is analytically given [3]. 3. Monte Carlo The kinetic approach outlined above neglects interac- tions fluctuations considering an (average) continuum loss of energy suffered by particles. In the case of protons, this approximation has a limited effect on the flux computation only at the highest energies (E > 100 EeV) [1, 7, 8]. In order to evaluate the effects of fluctuations on the expected nuclei flux, we have built a computation scheme alternative to the kinetic one, which uses the MC technique to simulate nuclei interactions. First of all, let us remark that fluctuations could be relevant only in the case of nuclei photo-disintegration. This follows from the fact that the pair-production process involving nuclei can be considered as an interaction process of the inside nucleon, therefore fluctuations in proton pair-production are irrelevant [8], and the same holds for nuclei. The SimProp MC simulation scheme that we have developed [4] is mono-dimensional: it does not take into account spatial distributions tagging sources only through their distance from the observer (red-shift). The MC simulation propagates particles in steps of red-shift following the injected nucleus, the secondary nuclei and protons produced at each photo-disintegration interaction and calculates their losses up to the observer, placed at red shift zero. The nuclear model on which SimProp is based is the same as is used for the kinetic approach (see [3, 4, and references therein]). The stochastic nature of the nuclei photo-disintegration process is modeled through the survival probability of a nucleus of atomic mass number A and Lorentz factor Γ P(Γ,z) = exp ( − ∫ z∗ z 1 τA(Γ,z ′) ∣∣∣∣ dtdz′ ∣∣∣∣dz′ ) (7) where z and z∗ are the values of the redshift of the current step (from z∗ to z). The SimProp code is designed in such a way that any red-shift distribution of sources and any injec- tion spectrum can be simulated. This is achieved by drawing events from a flat distribution in the red-shift of the sources and of the logarithm of the injection energy. Once the event is recorded at z = 0 the ac- tual source/energy distribution is recovered through a proper weight attributed to the event [4]. We will now compare the spectra obtained using SimProp [4] with the spectra calculated solving the kinetic equation associated to the propagation of nu- clei [3]. To pursue this comparison, a pure iron injec- tion with a power law injection of the type ∝ E−γg with γ = 2.2 has been assumed. The sources have been assumed to be homogeneously distributed in the red-shift range 0 < z < 3. In Fig. 1 the fluxes expected at z = 0 are shown for iron and secondary nuclei produced in the photo-disintegration chain suf- fered by primary injected irons. The points refer to the SimProp results, while the continuous lines refer to the fluxes computed in the kinetic approach. Good agreement between the two schemes is clearly visible in Fig. 1. At the highest energies the path-length of iron nuclei is very short. Therefore, to achieve good sampling in the MC simulation, higher statistics is needed; this is the reason for the larger errors bars in the SimProp results at the highest energies. Let us conclude by discussing why it is useful to go beyond the kinetic approach. The kinetic approach has the important feature of being analytical: the fluxes are computed mathematically by solving a first principles equation [3]. This means that the flux of pri- maries and secondaries is expressed in terms of several integrals that can be computed numerically, once the injection spectrum and the sources distribution are specified. In particular, the flux of secondary nuclei and nucleons produced in the photo-disintegration 705 Roberto Aloisio Acta Polytechnica chain of primary A0 is obtained by the numerical computation of A0 nested integrals and this computa- tion should be repeated each time the hypothesis on sources (injection and distribution) is changed. This computation, while it is always feasible numerically, takes some time. However the time can be substan- tially reduced by using a MC computation scheme. This follows from the fact that, as discussed above, it is possible within the SimProp approach to simulate different source distributions and injection spectra without repeating the overall propagation of particles. In this sense, a faster computation scheme is provided by the MC approach presented here, which is the minimal stochastic extension of the kinetic approach. Acknowledgements I am grateful to all collaborators with whom the works presented here were developed: V. Berezinsky, A. Gazizov and S. Grigorieva, from the theory group of the Gran Sasso Laboratory, D. Boncioli, A. Grillo, S. Petrera and F. Salamida, Auger group of L’Aquila University and P. Blasi, from the High Energy Astrophysics group of the Arcetri Astrophysical Observatory. References [1] R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, S. Grigorieva and B. Hnatyk, Astrop. Phys. 27 (2007) 76 [2] R. Aloisio and D. Boncioli, Astropart. Phys. 35 (2011) 152 [3] R. Aloisio, V. Berezinsky and S. Grigorieva (2012a), Astrop. Phys. in press DOI: 10.1016/j.astropartphys.2012.07.010, Astrop. Phys. in press DOI: 10.1016/j.astropartphys.2012.06.003 [4] R. Aloisio, D. Boncioli, A. Grillo, S. Petrera and F. Salamida, arXiv:1204.2970, (2012b), JCAP in press [5] Auger collaboration, Phys. Letters B 685 (2010) 239; Phys. Rev. Lett. 104 (2010) 091101 [6] V. Berezinsky, S. Bulanov, V. Dogiel, V. Ginzburg and V. Ptuskin, Astrophysics of Cosmic Rays, North-Holland, 1990 [7] V. Berezinsky, A. Gazizov and S. Grigorieva, (2006a) Phys. Rev. D 74 (2006) 043005 [8] V. Berezinsky, A. Gazizov and M. Kachelriess, (2006b) Phys. Rev. Lett. 97 (2006) 23110 [9] K. Greisen, Phys. Rev. Lett. 16 (1966) 748 [10] HiRes collaboration, Phys. Rev. Lett. 100 (2008) 101101 [11] HiRes collaboration, Phys. Rev. Lett. 104 (2010) 161101 [12] F. W. Stecker, M. A. Malkan and S. T. Scully, ApJ 648 (2006) 774 [13] G.T. Zatsepin and V.A. Kuzmin, Pisma Zh. Experim. Theor. Phys. 4 (1966) 114 Discussion Carlo Gustavino — The difference between Auger and the other experiments can be due to the fact they are looking from different hemispheres? Roberto Aloisio — This is an hypothesis that was recently put forward. I personally do not believe in such explanation because of the simple reason that at energies around (2 ÷ 3) × 1019 eV, where already the difference between Auger and HiRes starts, the universe visible in UHECR has a huge scale of the order of Gpc. Therefore it is very unlikely to have differences between observations carried out from the southern and northern hemispheres. 706 Acta Polytechnica 53(Supplement):703–706, 2013 1 Introduction 2 Kinetic Equations 3 Monte Carlo Acknowledgements References