CRPropa: a public framework to propagate UHECRs in the universe

To answer the fundamental questions concerning the origin and nature of ultra-high energy cosmic rays (UHECRs), it is important to confront data with simulated astrophysical scenarios. These scenarios should include detailed information on particle interactions and astrophysical environments. To achieve this goal one should make use of computational tools to simulate the propagation of these particles. For this reason the CRPropa framework was developed. It allows the propagation of UHECRs with energies $\gtrsim$10$^{17}$ eV and secondary gamma rays and neutrinos. The newest version, CRPropa 3, reflects an efficient redesign of the code as well as several new features such as time dependent propagation in three dimensions, galactic magnetic field effects and improved treatment of interactions, among other enhancements.


Introduction
The origin, nature and mechanisms of acceleration of ultra-high energy cosmic rays (UHECRs) are unanswered issues in astroparticle physics. To interpret the available experimental data one needs to model the source properties such as spectral index, mass composition and maximum acceleration energy, the distribution of sources, and the intervening cosmic magnetic fields. Effects arising from the interaction of cosmic rays with the pervasive photon backgrounds can also be relevant when constructing realistic scenarios.
Three observables measured in cosmic ray experiments are the spectrum, arrival directions and mass composition, the last one indirectly inferred from other observables, based on hadronic interaction models. They are affected by the interaction of UHECRs with photon fields, matter, as well as extragalactic and galactic magnetic fields. Any scenario aiming to elucidate the fundamental questions regarding the origin and nature of ultrahigh energy radiation should explain these observables simultaneously, taking into account propagation effects such as energy losses and magnetic deflections. Therefore, the development of computational tools to simulate realistic scenarios that fit the data is necessary. For this reason the CRPropa code [1][2][3][4] was created.
This article is organized as follows: in section 2 we briefly review the photon backgrounds present in the universe and the interactions taking place at ultra-high energies; in section 3 we describe the structure of the code; section 4 contains short descriptions of the new features a e-mail: rafael.alves.batista@desy.de of CRPropa 3; in section 5 we present some applications; and in section 6 we present the concluding remarks and outlook.

Interactions and photon backgrounds
Ultra-high energy cosmic rays lose energy during their propagation to Earth mainly through four processes: pair production, pion production, photodisintegration (in the case of nuclei) and adiabatic expansion of the universe.
The universe is permeated by photons with different wavelengths, which compose the extragalactic background light (EBL). The two main photon fields that should be taken into account in the propagation of UHE-CRs in the universe are the cosmic microwave background (CMB) and the cosmic infrared background (CIB). The first can be analytically estimated, whereas the second has to be obtained from observations. A third astrophysical background is the Universal Radio Background (URB). It has considerable effects only at 10 22 eV, or in the development of electromagnetic cascades, which can be inhibited depending on the density of radio photons.
The redshift evolution of the number density of the cosmic microwave background is given by where is the energy of the background photon. For the CIB the number density is determined through observations. Its redshift evolution is not trivial, so that a complete modeling of this problem requires measurements of the density of photons at different redshifts.
arXiv:1411.2259v1 [astro-ph.IM] 9 Nov 2014 The Journal's name The energy loss length for a nucleus of atomic number Z and mass A through production of electron/positron pairs can be written as [2] where α ≈ 1/137 is the fine structure constant, σ T is the Thomson cross section, h the Planck constant, k B the Boltzmann constant, Γ the Lorentz factor and f (Γ) a function taken from ref. [5]. The threshold energy for this interaction is E thr ≈ 5(meV/ ) EeV. Photopion production occurs when an EBL photon is scattered by a nucleon. In the case of protons, the two main interaction channels are The energy threshold for this process is E thr ≈ 70(meV/ ) EeV. For a CMB photon with ≈ 0.6 meV, E thr ≈ 4 × 10 19 eV, which is the expected energy for the well-known Greisen-Zatsepin-Kuzmin (GZK) cutoff. The pions produced through the interaction between nucleons and EBL photons decay as follows: This process is extremely important for multimessenger studies due to the production of secondary gamma rays and neutrinos. In CRPropa photopion production interactions are treated using the SOPHIA code [6]. The mean free path for this interaction for nuclei can be written as a combination of the ones for protons and neutrons. The interaction of atomic nuclei with EBL photons causes these nuclei to split into parts, through a photodisintegration process. In CRPropa photonuclear cross sections are obtained from the TALYS code [7]. The mean free path for this process can be written in terms of the cross section σ as follows: where Γ is the Lorentz factor, the photon energy, σ( ) the cross section of the nucleus-photon interaction, and max the maximum energy of the background photon, which is ∼10 meV for the CMB and ∼100 eV for the CIB. Unstable nuclei produced during photopion production or photodisintegration can have short lifetimes compared to the propagation length, and suffer decays during their trajectory. Our treatment of decays encompasses all relevant processes, namely α, β + and β − decays, and proton and neutron drippings, with tabulated lifetimes from NuDat 1 .
The expansion of the universe itself is another source of energy loss. This adiabatic energy loss is given by  Figure 1. Energy loss lengths for different processes: photopion production (orange), electron pair production (green), photodisintegration (purple), adiabatic expansion of the universe (gray) and total (black). Solid lines are for iron nuclei, and dashed lines for protons. The photon backgrounds used here are the cosmic microwave background and infrared background from ref. [8].
where E 0 is the initial energy. The energy loss length for photopion production, photodisintegration, pair production and adiabatic expansion of the universe are summarized in figure 1, for the case of iron and proton primaries.

Code structure
CRPropa 3 [3,4] is a reformulation of the previous versions of the code [1,2]. It is written in C++ with Python bindings. It inherited all features from CRPropa 2 and added new functionalities. The modular structure of this new version allows the user to extend the code for other applications. The modularity comes from the construction of the code, which handles separately sources, observers, particles and each interaction. Individual independent modules alters the property of the 'Candidate' class, which stores all details of the particle propagation, such as position, energy, type of particle, etc. These properties are updated at each step of propagation until a breaking condition is met or detection happens. CRPropa 3 supports shared memory parallel processing using OpenMP 2 , which allows fast simulations spanning a wide range of parameters to compare different scenarios.

Four-dimensional propagation
The simulation of UHECR propagation can be done in a one-dimensional (1D) environment, which allows the incorporation of cosmological effects such as the redshift dependence of the photon backgrounds, energy losses due to the adiabatic expansion of the universe, and source evolution, or in a three-dimensional (3D) environment, which can be used if one is interested in arrival directions, in addition to the spectrum and mass composition. In this case, cosmological effects cannot be taken into account because the information regarding the effective trajectory length of the particles, and therefore the redshift, is not known beforehand. The inclusion of cosmological effects in threedimensional simulations can be done by tracking the particles not only in the three spatial coordinates, but also in time, within a four-dimensional (4D) approach. This feature is now included in CRPropa 3.

Galactic Magnetic Field
CRPropa 3 allows the propagation of UHECRs considering several models of the galactic magnetic field (GMF). Implemented models include bissymmetric and antisymmetric spirals, as well as a toroidal field. The recent model by Jansson & Farrar [9,10] is also implemented, including regular, striated, and random turbulent components.
Any 3D or 4D simulation of the extragalactic propagation of UHECRs can be corrected for the effects of the galactic magnetic field. This is done a posteriori through a lensing technique first implemented in the PARSEC code [11], and adapted for CRPropa. For each energy there is a lens (matrix) which maps the directions of cosmic rays arriving at the edge of the galaxy to a direction observed at Earth. In this case no energy losses are considered.

Extragalactic Magnetic Fields
Realistic scenarios of source distribution and magnetic fields are important to constrain models of UHECRs. In this context, cosmological magnetohydrodynamical (MHD) simulations of the local universe can provide information regarding the matter distribution and intervening magnetic fields. In previous versions of CRPropa it was possible to consider a simple turbulent magnetic field, as well as magnetic fields from MHD simulations through uniformly spaced grids. The usage of multiresolution grids in MHD simulations has increased recently, and these higher resolution non uniform grids could be useful for cosmic ray propagation. For this reason CR-Propa 3 interfaces with external codes to handle MHD simulations from the SPH (Smooth Particle Hydrodynamics) code Gadget [12] and AMR (Adaptative Mesh Refinement) code RAMSES [13].

UHE photon propagation with the EleCa code
Ultra-high energy photons can be primary particles, emitted by an astrophysical object, or secondaries, generated through the interaction of primary cosmic rays with photon backgrounds. The interaction of these photons with the EBL induces the development of an electromagnetic cascade, which can be propagated within CRPropa with the external codes DINT [14], which calculates the photon spectrum by solving kinetic equations, or EleCa [15], a Monte Carlo code for propagating UHE photons in the universe.
The resulting photon spectrum calculated by DINT ranges from 10 8 eV up to 10 23 eV, while EleCa is restricted to higher energies ( 10 16 eV). Because DINT is based on transport equations, it is more efficient for lower energies, allowing a faster calculation of the spectrum compared to the more precise particle-by-particle Monte Carlo approach of EleCa. For now the propagation of photons with EleCa is done only in one dimension, and the effects of magnetic fields on the cascades are taken into account using a small angle approximation.

Updated photodisintegration cross sections
As mentioned before, in CRPropa 3 the treatment of photodisintegration is done using tabulated values for the cross section, taken from the TALYS code [7]. The previous version, CRPropa 2, uses photodisintegration cross sections from TALYS 1.0, whereas CRPropa 3 includes the upto-date TALYS 1.6 3 cross sections for A≥12, which improves the binning of the data around resonances, extends the available energy range, and enhances some reactions, among other improvements.
The impact of the new photodisintegration cross sections provided by TALYS 1.6 on the photodisintegration rates are approximately the same as the ones from TALYS 1.0 for the lightest and heaviest nuclei, but they differ significantly for intermediate masses. In figure 2 the impact of these differences on the interaction rate of UHECRs with EBL photons are shown for the case of carbon-12 and oxygen-16.  Figure 2. Photodisintegration interaction rates for 12 6 C (top) and 16 8 O (bottom) for TALYS 1.6 (green) and TALYS 1.0 (orange). 3 For detailed information on the changes made between TALYS versions 1.0 and 1.6, refer to the documentation in the website www.talys.eu. The Journal's name

Cosmic Infrared Background models
CRPropa 3 provides several options of models of the CIB. The user can choose between a set of predefined CIB models, namely the ones by Kneiske et al. [8], Franceschini [16], Dole et al. [17], Kneiske & Dole lower limit [18] and Stecker [19].

Applications
First we illustrate an application of CRPropa within a multimessenger approach, propagating UHECRs and obtaining the photon and neutrino counterparts. We assume a uniform distribution of sources with spectrum dN/dE ∝ E −2 . The infrared background model assumed is the one by Kneiske et al. [8]. The secondary photons and neutrinos spectra produced from the interaction of protons with background photons were also computed. The spectra for cosmic rays and secondaries are shown in figure 3. In this figure data from the Pierre Auger Observatory [20] are displayed for the sake of comparison.  We now present the results of a three-dimensional simulation. We use the large scale matter distribution from Dolag et al. [21], constrained in such a way to reproduce the observational data from cosmological surveys. This box, henceforth called matter distribution grid, is a cube grid of approximately 132 Mpc with a spacing of ∼300 kpc. Instead of using the magnetic field distribution from this MHD simulation, we use the one from Miniati [22], also used in several subsequent works [23,24] for the propagation of UHECRs. This model has a higher magnetic field strength compared to Dolag et al.. We use the profile magnetic field-density distributions from the Miniati simulation to obtain a relation between magnetic field and density. We then create a modulation grid by replicating the matter distribution grid, replacing the density in each cell by the corresponding magnetic field strength from the profile. The modulation grid has 256 3 cells covering a volume of approximately (132 Mpc) 3 . It is used to modulate another 256 3 grid with volume ∼(13.2 Mpc) 3 containing a realization of a turbulent Kolmogorov field with coherence length 500 kpc, periodically repeated to cover the complete simulation volume.

energy [eV]
To illustrate the propagation including effects of magnetic fields and matter distribution we arbitrarily assume a scenario composed of four species of atomic nuclei, hydrogen, helium, nitrogen and iron, with fraction 1, 0.5, 0.25 and 0.125, respectively. The sources follow the large scale distribution from Dolag et al. [21], and emit particles with a differential spectrum proportional to E −1.8 . The maximum rigidity of the source is R max = 10 19.8 eV. The maximum trajectory length considered was 2 Gpc.
In figure 4 the cosmic ray spectrum above 10 18 eV is shown.  The mass composition of the observed particles is expected to be different from the injected one, due to photodisintegration. In terms of experimental observables, a change in composition implies a change in the depth of the shower maximum ( X max ). In figure 5 the values of X max and σ(X max ) for this scenario are shown. They were estimated using the parametrization from refs. [25,26], assuming the EPOS-LHC hadronic interaction model [27].

energy [eV]
The skymaps containing the arrival directions of the simulated events are shown in figure 6 considering only extragalactic deflections and including effects of the galactic magnetic field according to the model by Jansson & Farrar [9,10].

Summary and outlook
The comparison between experimental data and theoretical models is crucial to understand the origin and nature of the UHECRs. Therefore, the development of computational tools that allow the detailed simulation of astrophysical scenarios including all relevant particle physics and astrophysical ingredients are essential to build realistic scenarios.
We have presented CRPropa 3, a public framework to propagate UHECRs,and secondary gamma rays and neutrinos in the universe. The code allows parallel processing, python steering and the inclusion of custom modules. The main new features are the four-dimensional propagation,  Figure 5. Estimated X max and σ(X max ) for the simulated data (orange line). Black markers correspond to data measured by the Pierre Auger Observatory [20], and gray regions are the systematics. Dotted lines are the predicted composition for pure proton and pure iron scenarios, according to the EPOS-LHC hadronic interaction model. galactic magnetic fields, new extragalactic magnetic field techniques, improved interaction tables and new CIB models. We have presented two applications for illustration purposes, and compared the results with measurements from the Pierre Auger Observatory. CRPropa 3 is in the final stages of development, is publicly available, and will soon be released. More information can be found on https://crpropa.desy.de.