CRPropa 3 – a Public Astrophysical Simulation Framework for Propagating Extraterrestrial UltraHigh Energy Particles
Abstract
We present the simulation framework CRPropa version 3 designed for efficient development of astrophysical predictions for ultrahigh energy particles. Users can assemble modules of the most relevant propagation effects in galactic and extragalactic space, include their own physics modules with new features, and receive on output primary and secondary cosmic messengers including nuclei, neutrinos and photons. In extension to the propagation physics contained in a previous CRPropa version, the new version facilitates highperformance computing and comprises new physical features such as an interface for galactic propagation using lensing techniques, an improved photonuclear interaction calculation, and propagation in time dependent environments to take into account cosmic evolution effects in anisotropy studies and variable sources. First applications using highlighted features are presented as well.
a,f]Rafael Alves Batista, a]Andrej Dundovic, b]Martin Erdmann, c]KarlHeinz Kampert, b,1]Daniel Kuempel,^{1}^{1}footnotetext: Corresponding author. b]Gero Müller, a]Guenter Sigl, a,d]Arjen van Vliet, b]David Walz, b,c,e]Tobias Winchen Prepared for submission to JCAP
CRPropa 3 – a Public Astrophysical Simulation Framework for Propagating Extraterrestrial UltraHigh Energy Particles

University of Hamburg,
II Institut für Theoretische Physik Luruper Chaussee 149, 22761 Hamburg, Germany 
RWTH Aachen University,
III. Physikalisches Institut A OttoBlumenthalStr., 52056 Aachen, Germany 
University of Wuppertal,
Department of Physics, Gaußstr. 20, 42097 Wuppertal, Germany 
now at Radboud University Nijmegen,
Department of Astrophysics/IMAPP, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands 
now at Vrije Universiteit Brussel,
Astrophysical Institute, Pleinlaan 2, 1050 Brussels, Belgium 
now at University of Oxford,
Department of Physics  Astrophysics, DWB, Keble Road, OX1 3RH, Oxford, UK
Keywords: ultra high energy cosmic rays, magnetic fields, ultra high energy photons and neutrinos
Contents
1 Introduction
The origin of high energy cosmic rays is still an open research question of fundamental interest for the understanding of our universe. Multiple aspects of cosmic rays have been investigated experimentally, most notably their steeply falling energy spectrum with a cutoff above EeV (eV) [1]. Cosmic ray arrival directions appear to be rather isotropic with a few exceptions only. As examples we mention a dipole signal above EeV [2], and a hot spot observed in the Northern hemisphere for energies EeV [3]. Concerning cosmic ray composition, measurements of the rate and shape of the depth of cosmic ray induced air showers reveal a composition with a contribution of large nuclear masses above EeV [4].
In recent years, high energy neutrinos have also been observed with a flux well above expectations from atmospheric showers [5]. Such extraterrestrial neutrinos with energies in the PeV regime (eV) are distributed all over the sky and may provide directional information on hadronic acceleration sites. So far, no significant autocorrelations, or correlations with matter distributions [5], or with cosmic ray arrival directions [6] have been observed.
Combining these experimental observations with current knowledge on large scale structures, magnetic fields and cosmic background fields leads to the following understanding of high energy cosmic rays. The lack of significant correlations of arrival directions above EeV with the galactic plane together with the limited cosmic ray confinement owing to the size of our galaxy and its magnetic field probably mean that these cosmic rays are of extragalactic origin. The exact transition from galactic to extragalactic cosmic rays is, however, not well understood yet [7, 8]. Taking advantage of the isotropic arrival distributions, at least bounds on the density of sources were derived depending on the cosmic ray energy [9].
The observed mass composition leads to predictions for deflection of the cosmic rays in the galactic and extragalactic magnetic fields that are larger by far than expected for protons. These deflections can reach tens of degrees for iron nuclei above EeV [10], and may be one of the reasons for the rather isotropic high energy cosmic ray sky. This could also explain the lack of correlations between neutrinos and cosmic rays.
Furthermore, attempts to explain the measured energy spectrum and composition data within a onedimensional astrophysical model prediction were undertaken where a nuclear composition has been adapted together with the slope of the injection spectrum at the sources and their maximum acceleration energy. The best fits favor injection spectra that are considerably harder than what is expected from shock acceleration theory. In addition, the preferred values for the maximally injected energies are comparatively low so that the “cutoff” of the spectrum may not be dominated by interactions with background fields as predicted by the socalled GZK effect [11, 12] but rather be caused by the limited source energies [13].
Obviously, the origin of cosmic rays is not explained easily but requires multiple aspects to be taken into account. Knowledge on many of these aspects has been acquired in the past decade or before, such that combining this knowledge in a numerical tool appears mandatory. The tool can then be used to develop different astrophysical predictions to be compared with various data distributions of different cosmic messengers. Such scenarios include models for the large scale distribution of the sources, their injection spectra and compositions, as well as models for the galactic magnetic field, and for the much less well known extragalactic magnetic fields and its structure. In this way experimental measurements can be maximally exploited to reject invalid scenarios, and to make progress in identifying an astrophysical scenario that is compatible with all measured distributions simultaneously.
The physics of nuclear decay and particle interactions with background fields is well known from laboratory experiments. Information on relevant background fields such as the cosmic microwave background (CMB) and the ultraviolet, optical, and infrared backgrounds (IRB) exists from astronomical observations, see e.g. [14]. Charged particle deflection in magnetic fields is precisely understood from electrodynamics. Progress in the description of the galactic magnetic field has been achieved recently by parametrizations respecting numerous Faraday rotation and starlight polarization measurements [16, 10, 15]. For extragalactic fields one can get information at least from simulations of the cosmological large scale galaxy structure which generate magnetic fields based on models of magnetohydrodynamics [17, 18]. Injection spectra at sources are not yet established, but can be obtained, e.g., from shock acceleration theory.
Furthermore, at energies around EeV and below, cosmic rays propagate over cosmological distances where the cosmological expansion, variability of low energy target photon backgrounds relevant for interactions, and deflection, even diffusion, in the structured cosmic magnetic fields all become relevant. Here information results from cosmological interpretations of astronomical observations [19].
The physics of secondary messengers, namely neutrinos and photons, is also well understood. Being weakly interacting matter particles, neutrinos have cross sections small enough to be essentially unaffected by background fields and magnetic fields. High energy photon interactions are precisely described by Quantum Electrodynamics and may produce lepton pairs or even electromagnetic showers in background fields.
Within this context, we have developed the program CRPropa version 3 as a general and versatile simulation tool that aims at efficient development of astrophysical predictions, and produces as output primary and secondary cosmic messengers such as protons, pions, nuclei, charged leptons, neutrinos, and photons. The program is publicly available. Its modular structure allows different components of a given astrophysical scenario to be combined and assembled. Users can also extend scenarios by including their own physics modules. In comparison to the previous version CRPropa 2 [20] most of the propagation physics implemented in CRPropa 3 is identical. However, the architecture and the code implementation have been completely reworked in order to optimally profit from modern programming design and computing techniques.
CRPropa 3 also contains new features, most notably models for the cosmological evolution of the infrared and radio backgrounds, corrections for taking into account the effects of cosmological expansion in the threedimensional modus that is used when simulating deflections in cosmic magnetic fields, evaluation of deflections in the galactic field, updates of the implemented photodisintegration processes, and a “fourdimensional mode” which allows to simulate time dependent scenarios by registering particle detections within a chosen time window. Furthermore, the development and propagation of electromagnetic cascades can now be simulated numerically using a Monte Carlo approach.
The paper is structured as follows. In section 2 we will briefly recap the inherited functionality and physics from the previous version CRPropa 2. Section 3 explains the novel code structure of the program before introducing the main new features for galactic and extragalactic propagation. The capabilities for simulating ultrahigh energy cosmic ray propagation in CRPropa 3 are demonstrated in section 4 in a few example applications. Finally, results are summarized in section 5.
2 Inherited Features from CRPropa 2
In the following we will introduce the features of CRPropa version 2 [20] that have been inherited by CRPropa 3. The publicly available software tool CRPropa 2 simulates the extragalactic propagation of UHE protons, neutrons, heavier nuclei and their secondary photons, electronpositron pairs and neutrinos. Included interactions are pair production, photopion production and photodisintegration with the cosmic microwave background (CMB) and the UV/optical/IR background (IRB) as well as nuclear decay. It can be run in onedimensional (1D) mode and in threedimensional (3D) mode. In 3D mode a 3D source distribution can be specified and deflections in extragalactic magnetic fields can be simulated as well. In 1D mode the only influence of the magnetic field that is taken into account is energy loss of pairs due to synchrotron radiation. The effects of cosmological, source and background light evolution with redshift can be included in the 1D mode as well. For the propagation of secondary photons and pairs a module called DINT is provided that solves the 1D transport equations for electromagnetic cascades initiated by electrons, positrons or gamma rays [21]. All these features of CRPropa 2 are available in CRPropa 3 as well.
3 New Features in CRPropa 3
3.1 Code structure and steering
An important step has been made in redesigning CRPropa 3 following actual standards for modular codes. Different aspects of the simulation (e.g. photonuclear interactions, boundary conditions, observer coordinates etc.) are separated into modules. Each module is independent of other modules and the probability, e.g. for an interaction, is calculated in each propagation step. To assure that the stepsize is small enough to process different modules in one propagation step in an arbitrary order, it can be limited by any of the modules. This default accuracy can be modified by the user if necessary. Since there are no direct dependencies between modules, any combination of modules can in principle be selected, allowing for multiple use cases and to study in detail individual propagation aspects. Therefore, each module can be replaced or added, making CRPropa 3 a flexible framework that can be extended without the need to modify other components. In this context, the simulation is implemented as an userdefined sequence of independent modules, which in turn modify objects of the main cosmic ray candidate class.
The cosmic ray candidate class incorporates the relevant information about the particle propagation, for instance the actual particle type, its comoving coordinates and velocities at different times and the list of secondary particles. The candidate properties are updated at each step of its propagation until a userprovided breaking condition is satisfied. Cosmic ray candidates can be created individually by the user or by means of a source model class, which takes the pertinent source properties, e.g. positions, spectrum, composition, and timedependencies as input.
A graphical visualization of the propagation process is given in figure 1. In this configuration, first the cosmic ray candidate is created by the source class. Then the modules sequentially process the cosmic ray candidate until the propagation is stopped by a user defined breaking condition. Output modules store all relevant information.
Cosmic ray propagation is a parallel task since interactions between cosmic rays are negligible. Current multicore processors can therefore be adequately utilized by just running multiple simulation instances in parallel. However, CRPropa 3 enables sharedmemory multiprocessing using OpenMP^{1}^{1}1http://openmp.org for easy multicore simulations. In this configuration, the ideal linear scaling is achieved up to about 8 threads. This limit is determined by critical sections in output modules and some external libraries (e.g. SOPHIA [22]).
CRPropa 3 is written in C++ and interfaced to Python using SWIG^{2}^{2}2http://www.swig.org. This allows the user to set up and steer simulations in a high level scripting language while all computations are performed with the underlying C++ code. The SWIG interface enables crosslanguage polymorphism, which can be used to extend a CRPropa simulation directly from the Python script that runs it. The user can for example write a custom simulation module in Python to be used in combination with the existing C++ modules. While the Python usage is the recommended steering mode, backward compatibility to the XML steering of CRPropa 2 is provided as well.
The performance of the code highly depends on simulation settings. For the simulation to finish, all candidates need to fulfill one of the given breaking conditions while iterating through all included modules in every propagation step. Hence, in 1D mode running time is determined solely by the execution time of the included modules, in 3D mode traversing magnetic fields dominates the overall performance: the larger deflections of cosmic rays in magnetic fields are, the smaller are the propagational steps resulting in more iterations of all modules per physical length and requiring more CPU time for the candidate to reach one of the spatial breaking conditions. To measure the performance and to find possible bottlenecks, one can use the PerformanceModule which calculates the percentage of time used in the total running time for every single module. More information about the performance in a specific scenario can be found in section 4.
3.2 Magnetic fields
As part of the restructuring, CRPropa 3 now supports any kind of magnetic field. The only requirement is the implementation of a getField function in the C++ MagneticField interface. This allows analytical fields, grid like fields, or fields from complex structures of any scale to be included. Table 1 summarizes the generic general and galactic magnetic field types implemented in CRPropa 3.
In addition, CRPropa 3 utilizes third party libraries to access data from adaptive mesh refinement and smooth particle simulations: The SAGA (SQLite AMR Grid Application) code^{3}^{3}3https://github.com/rafaelab/saga uses RTrees for fast access to RAMSES [26, 27] magnetohydrodynamical simulations data. Quimby^{4}^{4}4https://forge.physik.rwthaachen.de/projects/quimby [28] is a multiresolution grid library for fast access to huge compressed magnetic field grids and provides access to smooth particles from the GADGET [29, 30] file format.
The new structure also allows the implementation of magnetic field decorators, which can modify the given magnetic field on the fly. A periodic decorator turns any magnetic field into a periodic one, with the option of reflective boundaries. An evolution decorator allows for a cosmological scaling of type , where is the comoving magnetic field and the magnetohydrodynamic amplification of the field, cf. section 3.5. Multiple magnetic fields can also be grouped together and their magnetic fields can be superpositioned (added) in a list.
3.3 Galactic propagation
The propagation framework described above can, in principle, model cosmic ray propagation on any scale. When modeling the propagation of extragalactic cosmic rays entering the Milky Way with CRPropa 3 interactions can be neglected and only deflections in the magnetic field need to be considered. Instead of loading tabulated magnetic field data using the existing modules, dedicated modules describing the deflection in specific magnetic field models are used. The galactic magnetic field models proposed by Jansson & Farrar [10, 15] and Pshirkov et al. [16] are implemented and can be used as examples for other models. Forward and backward tracking of particles can be achieved by injection of the corresponding particles into the simulation.
However, to account for galactic magnetic field effects over large distances using forward propagation is computationally inefficient as most of the simulated particles miss the observer. As alternative, backtracking of cosmic rays with opposite charge from the observer to the edge of the galaxy is a much more efficient option to study possible trajectories of cosmic rays inside the Milky Way. In CRPropa 3 we provide an interface to the ‘lensing technique’ developed for the PARSEC software [31] that allows an efficient usage of backtracking simulations to account for the effects of deflection in the galactic magnetic field in forward simulations of extragalactic cosmic rays.
In the lensing approach, the trajectories of backtracked particles with rigidity are stored in matrices which are used to map discrete directions (pixel) indexed with outside the galaxy to discrete observed directions indexed with on Earth. The set of matrices form the ‘galactic lens’ that completely describes the deflection of extragalactic cosmic rays in the galactic magnetic field model. Using the matrices, a vector of the probabilities to observe a cosmic ray at energy from direction at the edge of the galaxy can be transformed by a matrix vector multiplication
(3.1) 
to obtain the probability distribution for energy . To avoid artificial distortions of the energy spectrum in this approach, all matrices have to be normalized by the maximum of unity norms of all matrices in the set.
The observed probability distributions can either be analyzed directly or be used to generate sets of individual simulated cosmic rays. Convenient interfaces to create the probability distributions from extragalactic CRPropa simulations and to sample data from the probability distributions are included in CRPropa 3.
The observed and injected directions are discretized using the HEALPix scheme [32]. The matrices are accessed in sparse compressed column major format using the Eigen^{5}^{5}5http://eigen.tuxfamily.org template library for linear algebra to minimize memory consumption and disk storage space while maximizing the computational performance for the matricesvector multiplications.
3.4 Photonuclear interactions
Photonuclear interactions with the CMB and IRB are the dominant energy loss processes for ultrahigh energy protons and nuclei. While the spectral shape of the CMB is well known, various models exist for the IRB. In addition to the Kneiske 2004 [33] model that was considered in CRPropa 2, the following IRB models are now available as well: Stecker 2005 [34], Franceschini 2008 [35], Finke 2010 [36], Dominguez 2011 [37] and Gilmore 2012 [14].
Furthermore, the implementation of photonuclear interactions has been improved in several ways. The pion production cross sections are now considered over a wider range of photon energies, and the integration procedure for calculating the interaction rates as well as the interpolation of the tabulated interaction rates were enhanced. This addresses the suggestions by Kalashev and Kido [38]. Additionally, the following improvements regarding photodisintegration were implemented. In CRPropa 2 photodisintegration cross sections were mainly obtained with TALYS 1.0 [39]; see ref. [20] for a detailed description. These cross sections are updated using the more recent TALYS version 1.6 (cf. ref. [40] for a full description of the changes in TALYS). More importantly, the parameters used in TALYS to model the giant dipole resonance have been adjusted to match the results given in [41]. The complete list of used parameters is given in appendix A, table 2. We verified that the recently published TALYS version 1.8 gives the same results with the used settings.
In CRPropa 3 the photonuclear cross sections are tabulated for photon energies in the rest frame of the nucleus of MeV in logarithmic steps of . This is done for all 169 isotopes in the range , with a lifetime of s. The branching ratios are taken into account for every channel that is computed in TALYS, namely with a simultaneous emission of up to eight nucleons in form of protons, neutrons, , , He and He nuclei. In practice, a large fraction of the resulting more than 25000 channels is of negligible impact for cosmic ray propagation. Thus, to increase performance, channels with branching ratios of less than 5% at every energy in the tabulated range are neglected, and the branching ratios of the remaining channels are scaled up accordingly. For the total cross section, however, all channels are considered. In contrast to the thinning procedure adopted in CRPropa 2, this procedure prevents a small systematic overestimation of the mean free path. For isotopes with mass numbers , the same set of cross sections as in CRPropa 2 is used, with the exception of Li, which is now modeled using the parametrization from Kossov [42]; as well as omitting He and Li, which, due to their short lifetime , exhibit negligible photodisintegration. As a result, TALYS is not used for any isotope with , as is recommended in ref. [40]. In total, photodisintegration is considered for 183 isotopes and 2200 channels.
Alternatively, CRPropa 3 provides the option of modeling photodisintegration for all isotopes using the parametrization from Kossov [42]. While the Kossov parametrization models both photodisintegration and photohadronic interactions, only the photodisintegration part up to the pion production threshold is considered here, since pion production is treated separately in CRPropa. Also, since the parametrization only models the total cross section, branching ratios are computed with TALYS as described above. A comparison of the two approaches to the available measured nuclear cross sections from the IAEA handbook on photonuclear data [43] indicates a slightly better agreement with the TALYS version. The differences in typical UHECR scenarios is at level of 10% in the spectral flux; see ref. [44] for more details.
3.5 Cosmological effects
Cosmological effects can affect the propagation of UHECRs. The adiabatic expansion of the universe implies a reduction in the momentum of a cosmic ray by a factor , where is the redshift whose evolution is given by
(3.2) 
in the standard CDM cosmology. Here, the Hubble parameter at present time is , is the density of matter in the universe, encompassing both baryonic and dark matter, and is the cosmological constant, assuming a flat universe () [19].
Interaction rates for different processes depend on the redshift at which the cosmic ray interacts with the photon backgrounds, because the number density and spectral shape of these radiation fields are evolving with time. In the case of the CMB, the spectral number density evolves passively as . On the other hand, the IRB is determined by the sum of galactic luminosities integrated over the entire age of the universe, hence its dependence is nontrivial. Similarly to CRPropa 2, we consider the evolution of the IRB through a global scaling factor , thereby assuming the spectral shape of the number density to be constant. The scaling factor is obtained for each IRB model (cf. section 3.4) as the integral over the comoving spectral number density, relative to its value at the present time. Using this approximation, the interaction rate is given by
(3.3)  
(3.4) 
In contrast to CRPropa 2, the integral for the IRB is performed over the entire modeled energy range, instead of just the part where the IRB dominates over the CMB. We have also tested a more detailed treatment, in which the exact redshift evolution of the IRB is taken into account. We found that this affects the propagated spectra due to the photopion production by less than 1% in realistic source scenarios such as the one shown in sec. 4.2. In the more extreme case of high redshift sources () and only photopion production in the IRB, this value is below 10%. This uncertainty for photodisintegration is expected to be roughly the same, although this case was not explicitly tested.
Including cosmological effects in the simulation requires an a priori knowledge of the propagation length, thus the redshift at the time of emission. In a onedimensional environment this is trivial, since the redshift corresponds to the distance of the source to the observer. In a threedimensional environment, however, the situation is more complicated because the effective propagation length can change due to deflections caused by intervening magnetic fields, and due to the redshift dependence of the photon backgrounds. To take into account simultaneously cosmological effects as well as deflections by magnetic fields, in CRPropa 3 a generalization of the 3D tracking of particles encompassing the time, respectively redshift, dimension is introduced, which is henceforth called ‘4D mode’. This new feature extends the notion of a threedimensional observer to four dimensions, so that the observer is considered a hypervolume composed by a sphere of a given radius and a redshift window of size around . The calculations for is obtained by extrapolating the corresponding quantities down to the desired value.
The redshift window (, corresponding to ) has to be smaller than all relevant time scales on which the cosmic flux could vary, unless one is interested in studies of time variability of sources, in which case the window should be as small as possible. This means that if the Hubble time is expressed by , then is already a sufficient constraint if the density of background photons vary slowly with redshift. Moreover, at the extreme high energy part of the spectrum () we are limited by interaction horizons, which are smaller than , corresponding to a small redshift (). For adiabatic losses are much larger than the others, and that is the region where the effects of the redshift window would matter. We have performed tests to check the convergence of the spectrum for different in a scenario with null magnetic field. The differences between and are of the order of 10 at 1 EeV and 10 at 0.1 EeV.
CRPropa 3 uses comoving coordinates internally. This implies that the flux dilution of comoving magnetic fields with redshift is implicitly taken into account. Additionally, a redshift scaling of the form can be applied to any magnetic field to account for a general field damping resulting in a total magnetic field evolution of . An explicit evolution of structured extragalactic magnetic field models resulting from magnetohydrodynamical simulations is currently not implemented, but can be added as described in section 3.2.
The use of comoving coordinates also implies that the comoving density of an arbitrary distribution of sources is kept constant. The redshift at the time of emission can be set by hand, or randomly picked according to a source evolution of type , where e.g. the star formation rate for is typically described with [45]. Note, that the evolution parameter for the source density and magnetic field need not be the same, and are completely unrelated.
3.6 Secondary messengers
Neutrinos and gamma rays produced in UHECR propagation are important messengers as well. Neutrinos result from the decay of the charged pions produced in photopion production and in nuclear decays. Being weakly interacting particles they are only subject to adiabatic energy loss and propagate on straight lines. Compared to the previous versions of CRPropa neutrinos are now processed directly by the module chain, allowing an analysis of secondary neutrinos by using the same observer objects as for the UHECR nuclei in a single simulation chain.
High energy gamma rays are produced in the decay of neutral pions and by inverse Compton scattering of high energy electrons with background photons. Secondary photons from photonuclear interactions that result from secondary nuclei decaying from excites states are not considered for now [46]. Electrons are created in electronpositron production, beta decay of nuclei, and the decay of muons from the decay of charged pions. Together they form an electromagnetic cascade, in which the high energy photons dominantly interact with background photons via pair production and double pair production . The electrons lose energy via inverse Compton scattering , triplet pair production and synchrotron radiation from gyrating in magnetic fields.
To calculate the development of the electromagnetic cascade CRPropa 3 provides interfaces to shipped versions of the specialized codes DINT [21], which was also part of the previous versions of CRPropa, and as a new feature also to EleCa [47]. DINT calculates the observed spectra based on the transport equations given in ref. [21]. It is therefore computationally efficient and thus particularly suited for calculations at lower energies where the cascade consists of many particles. Conversely, EleCa provides a full Monte Carlo tracking of the individual particles with stochastic treatment of the energy losses. Both codes currently focus on 1D propagation allowing the investigation of the resulting energy spectra of secondary particles.
The shipped version of the EleCa code has been improved regarding its performance resulting in a speed up of a factor of 3.6 compared to the baseline in a benchmark application. The speed up was achieved by replacing the onthefly calculation of the differential cross sections for the individual processes by interpolation between precalculated values and by code optimizations. The relative difference in the differential cross section between the onthefly calculation and the interpolation is smaller than and thus negligible.
To access the cascade codes, the secondary photons and electrons generated in ultrahigh energy nuclei propagation are written to a separate file which can then be further processed outside of the module chain by DINT or EleCa. Additionally, we provide an interface to a combined propagation in which particles above a customizable threshold energy are propagated individually with EleCa, while particles below the threshold energy are processed with DINT, cf. section 4.1.
4 Example Applications Focussing on New Features
4.1 1D simulation including secondary particles
To demonstrate the new features of CRPropa 3 we investigate the production of secondary messengers generated in the propagation of UHECR proton and iron primaries in a 1D simulation. We inject 10000 primaries from homogeneously distributed sources from a minimum distance of 3 Mpc up to a maximum distance of 1000 Mpc. The energies of the primaries are selected with a spectrum following a power law with index between a minimum energy of 1 EeV and a maximum energy of 1000 EeV. We include all available energy losses as described in the previous section using the default photon background models in the simulation.
The energy distributions of the secondaries injected into the electromagnetic cascade are shown in figure 2 separated into generating processes and type of primary.
For protons, the majority of secondary particles are low energy electrons created in electron pair production. Higher energy electrons are also generated in the decay of neutrons, which were created in photomeson production. The highest energetic secondaries are photons from the decay of the neutral pions. For iron primaries, the created secondaries are also electrons from pair production and high energy photons from the decay of neutral pions. However, the number of photons is lower and the photons have a lower mean energy due to the lower energy per nucleon of the particles. The number of created electrons from pair production, however, is higher than in the proton case due to the larger number of secondary nuclei that are created from photodisintegration. No secondaries from decay of radioactive nuclei are generated in this simulation.
The secondary photons and electrons are then propagated with DINT and the combined DINTEleCa propagation module using a threshold energy of 0.5 EeV. The observed spectra are shown in figure 3 as histograms for the combined propagation and gray surface for the DINT only propagation.
For energies above a few hundred TeV (eV) the results of both propagation modules are in reasonable agreement for iron primaries, considering the relative low number of high energy photons in this simulation. For energies below a few hundred TeV the propagation with EleCa for the highest energetic photons and electrons results in a decreased spectrum for proton and iron primaries compared to the DINT only propagation. Furthermore, in case of proton primaries also differences between both modules above a few hundred PeV are visible.
The inclusion of secondary electrons increases the size of the simulation output significantly, amounting to 2.5 GB of uncompressed output in ascii format in case of the proton primaries and 23 GB for the iron nuclei. The secondary electrons are important for the estimation of the photon flux below approximately 5 EeV.
To benefit from the new modular and flexible structure we are currently developing a photon and electron propagation module for simulating electromagnetic cascades which will be integral part of the CRPropa 3 framework. The results will be reported separately.
The resulting flux of secondary neutrinos is displayed in figure 4. In the case of proton primaries, approximately four times more neutrinos are generated in decays of compared to decays. In the case of iron primaries, approximately ten times more neutrinos are generated from decay than from the decay of . For the same number of injected particles, approximately twice the number of neutrinos are observed for iron primaries than for proton primaries. However, in this case twenty times more UHECR are observed from the same number of injected primaries.
4.2 3D simulation including extragalactic and galactic deflections
To show the capabilities of the 3D mode of CRPropa 3 an example is presented here including deflections in extragalactic and galactic magnetic fields as well as a source distribution following the largescale matter density. In this scenario the sources are distributed randomly with a source density of following the largescale structure formation simulations from Dolag et al. [48], which have been constrained by the observed largescale density field with a radius of 110 Mpc around the Milky Way.
The structure of the implemented extragalactic magnetic field model in this example is obtained from the same constrained simulations by Dolag et al. The overall magnetic field strength, however, is smaller in their simulations than in the similar but unconstrained simulations done by Sigl et al. [49]. Sigl et al. scaled the magnetic field strength such that the magnetic field in the core region of a Comalike galaxy cluster is of the same order as indicated by Faraday rotation measures. To emphasize cosmicray deflections, we derived the magnetic field strength for this example from the relation between matter density and magnetic field strength obtained from the simulations by Sigl et al. This has been done by assigning to each cell in the largescale matter density grid from Dolag et al. the magnetic field strength from the simulations by Sigl et al. corresponding to the matter density in that cell. One thus obtains a magnetic field modulation grid (see also section 3.2) of cells covering a volume of . A higherresolution magnetic field with Fourier modes taken from a Gaussian distribution with given by a Kolomogov power spectrum and random polarization, with a coherence length of 500 kpc and a total volume of , is then periodically repeated to cover the complete modulation grid. Cosmic rays with trajectory lengths up to 4 Gpc are taken into account by employing reflective boundary conditions when they reach the edge of the simulation box.
Apart from magnetic field deflections, all available photonuclear and decay interactions (see section 3.4) on both the CMB and the IRB have been included in the extragalactic propagation as well. In this 3D simulation the photon background is taken as timeindependent and adiabatic losses are neglected. The implemented IRB model is that from Gilmore 2012 [14]. Two different scenarios are presented here, one for pure proton injection at the sources and one for pure iron injection. In both cases the cosmic rays are initiated at their sources following a power law spectrum with a broken exponential cutoff:
(4.1) 
with the number of injected particles and the particle energy. For this example the particles are injected with a spectral index of and cutoff energy EeV down to a minimum energy of EeV.
We consider the observer to be a sphere with a given radius , which should be small enough to avoid spurious effects arising from fluctuations in the magnetic fields in different regions of the observer sphere. The smaller is, the closer the simulated observer is to reality, but the larger the required computational time for a sufficient number of observed events. A good tradeoff is to take , where is the typical scattering length for a charged particle in the magnetic field. This condition ensures that no significant diffusion will occur within the observer. From numerical studies of the trajectories of protons in this magnetic field we have obtained . For this reason we set .
Due to the finite size of the observer, another important factor that should be taken into account is the multiplicity of detections. In CRPropa it is possible to consider multiple detections of the same particle. However, in reality this would not happen and simulating multiple detections of the same cosmic ray is not in accordance with Liouville’s theorem. To solve this issue we randomly select one hit among all detections of the same particle. Therefore, when multiple hits within the same periodic box from the same original cosmic ray are registered, one of these hits is chosen randomly and only that hit is used in the analysis.
For the galactic propagation the lensing technique described in section 3.3 has been used. The implemented galactic magnetic field model is the Jansson and Farrar 2012 model [10, 15] including the random largescale and turbulent smallscale component. Resulting sky maps before and after galactic magnetic field deflections for proton and iron injections at the sources are shown in figure 5. These sky maps have been normalized so that the bin with the maximum number of hits is set to one. The range of the color scale has been set to values from 0.3 to 0.7 to better visualize the differences between the sky maps.
The coordinateindependent angular power spectra,
(4.2) 
with multipoles for these sky maps are shown in figure 6, normalized to . The sky maps and power spectra given here are for the full energy range ( EeV) and, as the flux is decreasing with energy, are dominated by cosmic rays with energies EeV. In this energy range the iron injection scenario consists mainly of light nuclei, products from photodisintegrated iron nuclei, and is therefore similar to the proton injection scenario. The dipole amplitude (see the first bin of figure 6 for ) for the proton injection scenario before (after) deflections in the GMF is (0.066) and for the iron injection scenario is (0.069). For comparison the expected confidence level upper bounds that would result from fluctuations of isotropic distributions with , and events are shown. This indicates the sensitivity of a hypothetical experiment with full sky coverage to detect the level of anisotropy presented in this example application for different numbers of detected events with energies above 1 EeV. The energy spectra and compositions at the observer are depicted in figures 7 and 8 together with the spectra and compositions for the 4D cases (cf. section 4.3).
4.3 Cosmological effects using 4D simulations
As an application of the redshift dependent 4D mode we use a similar setup to the 3D example, but now including cosmological effects. Adiabatic energy losses due to the expansion of the universe are included. Events are detected when they reach the observer within a redshift window of size , i.e. only events that are recorded in the redshift range are taken into account. For efficiency purposes we choose an observer of radius 1 Mpc, unlike in the 3D example, in which a 100 kpc observer is used. This does not significantly affect the spectrum and composition because the difference between the observer sizes in these two cases is small.
We adopt the same settings as in the 3D example and simulate two scenarios: pure iron and pure proton compositions. The maximum energy and source spectral index are the same as before. The only difference is the additional energy loss process (adiabatic losses) and the redshift dependence of the CMB and IRB. The sources in the 4D scenario are assumed to have a uniform redshift distribution up to . The spectra for both scenarios are shown in figure 7, and the mean and variance of the logarithm of the mass composition, and , observed at Earth for the iron cases are shown in figure 8, together with the universal spectrum, i.e., the spectrum obtained from a onedimensional simulation assuming a uniform distribution of sources.
In figure 7 we notice that the differences in the spectrum are overall small, which suggests that the analyzed scenarios are fairly close to the universal case. Nevertheless, there are noticeable effects in the observed composition for the iron source scenario as illustrated in figure 8. Note that the common feature of an abrupt increase in average mass at around is due to kinematics: An iron nucleus injected with an energy of several hundred EeV typically suffers complete photodisintegration into 56 nucleons of energy each. Due to the hard source spectrum these secondary protons dominate the observed flux up to the energy that corresponds to the cutoff in the source spectrum of the initial iron nuclei. One can see that for in the 3D case is smaller than in the 4D case. This is explained by the effect of adiabatic energy losses and generally stronger interactions in the 4D case, which are neglected in the 3D simulation. Indeed, additional energy losses at highest energies lower the energy up to which secondary protons can contribute and the transition is less abrupt. Similar reasoning holds for .
In the absence of intervening magnetic fields the composition as well as the spectrum for the 1D scenario are expected to be the same as for 4D. In this example, however, for , the 1D scenario has a slightly larger with respect to 4D. This can be explained by the presence of intervening magnetic fields in combination with a discrete source distribution, which cause an increase in the effective propagation time of the particles, and hence the probability of an interaction to occur during propagation. Another effect that, although subdominant in this example, might play a role at low energies () is the magnetic horizon, which spawns a suppression in the flux of cosmic rays when the trajectory lengths become comparable to the age of the universe [50, 51, 52, 53].
Although we have taken redshift effects into account, the comoving magnetic field can also evolve with redshift approximately as , as discussed in section 3.5. We have not considered the redshift evolution of the field, which is equivalent to taking in the aforementioned equation. The parameter is expected to be nonzero if the magnetic field is significantly affected by magnetohydrodynamic processes taking place during structure formation. Moreover, this approximation is adequate for very high energies (), where cosmic rays most likely originate in the nearby universe. At energies of the order of a few EeV the redshift evolution of the magnetic field might play an important role. We note, however, that even though the redshift evolution of the magnetic field can be taken into account in CRPropa this is outside the scope of the present work.
5 Summary
The simulation of galactic and extragalactic cosmic ray propagation plays an essential role in understanding astrophysical processes at ultrahigh energies. In this paper we introduced the new version of the publicly available cosmic ray propagation code CRPropa 3. To interpret the data collected by largescale cosmic ray observatories, the code was completely rewritten to be flexible enough to cover the large parameter space of astrophysical scenarios. The modular structure, included in version 3, enables to combine independent modules to study multiple use cases and even to extend the code by new individually specified modules.
While inheriting all features from the previous version, CRPropa 3 introduces additional functionalities. As a consequence of the modular and flexible code structure any kind of magnetic field is now supported and additionally modelization of deflection in galactic magnetic fields has been improved. Furthermore, the lensing technique provides a computationally efficient method to calculate trajectories of cosmic rays inside the Milky Way and can be applied for the most commonly used galactic magnetic field models.
The calculation of photonuclear interactions has been updated to latest models and a number of new IRB models have been implemented. In the case of 3D simulations cosmological effects are now taken into account and the evolution of the IRB with redshift is now implemented with a redshift dependent scaling. To take into account the a priori unknown propagation length and the resulting changes in photon background, CRPropa 3 extends the notion of a 3D observer to a 4D mode by taking the redshift (time) into account.
To complete the multimessenger picture, CRPropa 3 includes the production and propagation of secondary particles such as neutrinos and gamma rays. In addition to the transport code DINT for electromagnetic cascades we now incorporated a full Monte Carlo code to calculate the electromagnetic cascade at energies above 0.1 EeV based on the EleCa code. A combination of the two codes is available enabling to follow cascades down to TeV.
Information on downloading the code, usage and example applications can be found at http://github.com/CRPropa/CRPropa3. Questions and comments can be submitted to the ticketing system https://github.com/CRPropa/CRPropa3/issues. The features of CRPropa 3 can also be explored online without installation at https://vispa.physik.rwthaachen.de.
Appendix A Giant dipole resonance parameters
Photodisintegration is the most important interaction for cosmic ray nuclei with energies eV. Cross sections for this interaction are dominated by the giant dipole resonance for photons with energies MeV in the nucleus rest frame. In ref. [41] an “accurate description” of the available experimental data was found using a preliminary version of TALYS. TALYS was used in this comparison [54] with the giant dipole resonance parameters from the IAEA atlas [43]. In contrast, the publicly available versions of TALYS by default uses the giant dipole resonance parameters from the RIPL2 database [55], and the predicted cross sections are in poor agreement with the experimental data. Thus, in CRPropa 3 TALYS is used with the giant dipole resonance parameters of the IAEA atlas, if available, as the resulting cross sections are in much better agreement with the available measurements. In figure 9 we show the total photodisintegration cross sections for C and Si. The Kossov [42] parametrization and TALYS 1.6 with adjusted giant dipole resonance parameters are in reasonable agreement with experimental data. Both models are implemented in CRPropa 3. The complete list of used giant dipole resonance parameters is given in Table 2.
Isotope  [MeV]  [mb]  [MeV]  [MeV]  [mb]  [MeV]  Source 

C  22.70  21.36  6.00  Atlas  
N  22.50  27.00  7.00  Atlas  
O  22.35  30.91  6.00  Atlas  
Na  23.00  15.00  16.00  Atlas  
Mg  20.80  41.60  9.00  Atlas  
Al  21.10  12.50  6.10  29.50  6.70  8.70  RIPL2 
Si  20.24  58.73  5.00  Atlas  
Ar  20.90  50.00  10.00  Atlas  
Ca  19.77  97.06  5.00  Atlas  
V  17.93  53.30  3.62  20.95  40.70  7.15  RIPL2 
Mn  16.82  51.40  4.33  20.09  45.20  4.09  RIPL2 
Acknowledgments
The authors would like to thank Carmelo Evoli for contributions to the code, Lukas Merten for helpful comments and suggestions and Christopher Heiter and Mariangela Settimo for improvements in the photon propagation calculation. Thanks also to Marcus Wirtz for improving the galactic magnetic field modeling. We are grateful to Arian Koning and Stéphane Goriely for guidance in using the TALYS package. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the collaborative research centre SFB 676, by the Bundesministerium für Bildung und Forschung (BMBF) and by the Helmholtz Alliance for Astroparticle Physics (HAP) funded by the Initiative and Networking Fund of the Helmholtz Association including the cost takeover of the open access processing charge. AvV acknowledges financial support from the NWO astroparticle physics grant WARP.
References
 [1] Pierre Auger Collaboration, J. Abraham et al., Measurement of the energy spectrum of cosmic rays above eV using the Pierre Auger Observatory, Phys. Lett. B 685 (2010) 239 [arXiv:1002.1975].
 [2] Pierre Auger Collaboration, A. Aab et al., Large Scale Distribution of Ultra High Energy Cosmic Rays Detected at the Pierre Auger Observatory With Zenith Angles up to , Astrophys. J. 802 (2015) 111 [arXiv:1411.6953].
 [3] Telescope Array Collaboration, R. U. Abbasi et al., Indications of IntermediateScale Anisotropy of Cosmic Rays with Energy Greater Than 57 EeV in the Northern Sky Measured with the Surface Detector of the Telescope Array Experiment, Astrophys. J. 790 (2014) L21 [arXiv:1404.5890].
 [4] Pierre Auger Collaboration, A. Aab et al., Depth of maximum of airshower profiles at the Pierre Auger Observatory. II. Composition implications, Phys. Rev. D 90 (2014) 122006 [arXiv:1409.5083].
 [5] IceCube Collaboration, M. G. Aartsen et al., Observation of HighEnergy Astrophysical Neutrinos in Three Years of IceCube Data, Phys. Rev. Lett. 113 (2014) 101101 [arXiv:1405.5303].
 [6] IceCube and Pierre Auger and Telescope Array Collaborations, M. G. Aartsen et al., Search for correlations between the arrival directions of IceCube neutrino events and ultrahighenergy cosmic rays detected by the Pierre Auger Observatory and the Telescope Array, JCAP 001 (2016) 037 [arXiv:1511.09408].
 [7] G. Giacinti, M. Kachelriess, D. V. Semikoz and G. Sigl, Cosmic Ray Anisotropy as Signature for the Transition from Galactic to Extragalactic Cosmic Rays, JCAP 007 (2012) 031 [arXiv:1112.5599].
 [8] Pierre Auger Collaboration, P. Abreu et al., Constraints on the origin of cosmic rays above eV from large scale anisotropy searches in data of the Pierre Auger Observatory, Astrophys. J. 762 (2012) L13 [arXiv:1212.3083].
 [9] Pierre Auger Collaboration, P. Abreu et al., “Bounds on the density of sources of ultrahigh energy cosmic rays from the Pierre Auger Observatory,” JCAP 1305 (2013) 009 [arXiv:1305.1576].
 [10] R. Jansson and G. R. Farrar, A New Model of the Galactic Magnetic Field, Astrophys. J. 757 (2012) 14 [arXiv:1204.3662].
 [11] K. Greisen, End to the CosmicRay Spectrum?, Phys. Rev. Lett. 16 (1966) 748.
 [12] G. T. Zatsepin, & V. A. Kuz’min, Upper limit of the spectrum of cosmic rays PZETF 4 (1966) 114.
 [13] Pierre Auger Collaboration, A. di Matteo et al., Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory, Proc. 34th International Cosmic Ray Conference (ICRC), July 2015, The Hague, The Netherlands, [arXiv:1509.03732].
 [14] R. C. Gilmore, R. S. Somerville, J. R. Primack and A. Dominguez, Semianalytic modeling of the EBL and consequences for extragalactic gammaray spectra, Mon. Not. Roy. Astron. Soc. 422 (2012) 3189 [arXiv:1104.0671].
 [15] R. Jansson and G. R. Farrar, The Galactic Magnetic Field, Astrophys. J. 761 (2012) L11 [arXiv:1210.7820].
 [16] Pshirkov, M. S., P. G. Tinyakov, P. P. Kronberg and K. J. NewtonMcGee, Deriving the Global Structure of the Galactic Magnetic Field from Faraday Rotation Measures of Extragalactic Sources, Astrophys. J. 738 (2011) 192 [arXiv:1103.0814].
 [17] K. Dolag and F. A. Stasyszyn, An MHD Gadget for cosmological simulations, Mon. Not. Roy. Astron. Soc. 398 (2009)1678 [arXiv:0807.3553].
 [18] F. Miniati and D. F. Martin, Constrainedtransport Magnetohydrodynamics with Adaptive Mesh Refinement in CHARM, Astrophys. J. Suppl. 195 (2011) 5 [arXiv:1103.1878].
 [19] Particle Data Group, K. A. Olive et al., Review of Particle Physics, Chin. Phys. C 38 (2014) 090001.
 [20] K.H. Kampert et al., CRPropa 2.0  a public framework for propagating high energy nuclei, secondary gamma rays and neutrinos, Astropart. Phys. 42 (2013) 41 [arXiv:1206.3132].
 [21] S. Lee, On the propagation of extragalactic highenergy cosmic and gammarays, Phys. Rev. D 58 (1998) 043004 [arXiv:astroph/9604098].
 [22] A. Mucke, R. Engel, J. P. Rachen, R. J. Protheroe and T. Stanev, SOPHIA: Monte Carlo simulations of photohadronic processes in astrophysics, Comput. Phys. Commun. 124 (2000) 290 [astroph/9903478].
 [23] G. Giacinti, M. Kachelrieß, D. V. Semikoz, and G. Sigl, Ultrahigh energy nuclei in the turbulent Galactic magnetic field, Astropart. Phys. 35 (2011) 192 [arXiv:1104.1141].
 [24] M. Prouza and R. Smida, The galactic magnetic field and propagation of ultrahigh energy cosmic rays, Astron. Astrophys. 410 (2003) 1 [astroph/0307165].
 [25] X. H. Sun, W. Reich, A. Waelkens and T. Enßlin, Radio observational constraints on Galactic 3Demission models, Astron. Astrophys. 477 (2008) 573 [arXiv:0711.1572].
 [26] R. Teyssier, S. Fromang and E. Dormy, Kinematic dynamos using constrained transport with high order godunov schemes and adaptive mesh refinement, J. Comput. Phys. 218 (2006) 44 [astroph/0601715].
 [27] S. Fromang, P. Hennebelle and R. Teyssier, A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics,Astron. & Astrophys. 457 (2006) 371 [astroph/0607230].
 [28] G. Müller, Static Multiresolution Grids with Inline Hierarchy Information, [arXiv:1512.03172].
 [29] V. Springel, N. Yoshida and S. D. M. White, GADGET: A Code for collisionless and gasdynamical cosmological simulations, New Astron. 6 (2001) 79 [astroph/0003162].
 [30] V. Springel, “The Cosmological simulation code GADGET2,” Mon. Not. Roy. Astron. Soc. 364 (2005) 1105 [astroph/0505010].
 [31] H. P. Bretz, M. Erdmann, P. Schiffer, D. Walz and T. Winchen, PARSEC: A Parametrized Simulation Engine for UltraHigh Energy Cosmic Ray Protons, Astropart. Phys. 54 (2014) 110 [arXiv:1302.3761].
 [32] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke and M. Bartelman, HEALPix  A Framework for high resolution discretization, and fast analysis of data distributed on the sphere, Astrophys. J. 622 (2005) 759 [astroph/0409513].
 [33] T. M. Kneiske, T. Bretz, K. Mannheim and D. H. Hartmann, Implications of cosmological gammaray absorption. 2. Modification of gammaray spectra, Astron. Astrophys. 413 (2004) 807 [astroph/0309141].
 [34] F. W. Stecker, M. A. Malkan and S. T. Scully, Intergalactic photon spectra from the far IR to the UV lyman limit for and the optical depth of the universe to high energy gammarays, Astrophys. J. 648 (2006) 774 [astroph/0510449].
 [35] A. Franceschini, G. Rodighiero and M. Vaccari, Extragalactic opticalinfrared background radiation, its time evolution and the cosmic photonphoton opacity, Astron. & Astrophys. 487 (2008) 837 [arXiv:0805.1841].
 [36] J. D. Finke, S. Razzaque and C. D. Dermer, Modeling the Extragalactic Background Light from Stars and Dust, Astrophys. J. 712 (2010) 238 [arXiv:0905.1115].
 [37] A. Domínguez et al., Extragalactic background light inferred from AEGIS galaxySEDtype fractions, MNRAS 410 (2011) 2556 [arXiv:1007.1459].
 [38] O. E. Kalashev and E. Kido, Simulations of Ultra High Energy Cosmic Rays propagation, J. Exp. Theor. Phys. 120 (2015) 5 790 [arXiv:1406.0735].
 [39] S. Goriely, S. Hilaire and A. J. Koning, Improved predictions of nuclear reaction rates with the TALYS reaction code for astrophysical applications, Astron. & Astrophys. 487 (2008) 767 [arXiv:0806.2239].
 [40] A. J. Koning, S. Hilaire and M. Duijvestijn, TALYS1.6 user manual, www.talys.eu
 [41] E. Khan, S. Goriely, D. Allard, E. Parizot, T. Suomijarvi, A. J. Koning, S. Hilaire and M. C. Duijvestijn, Photodisintegration of ultrahighenergy cosmic rays revisited, Astropart. Phys. 23 (2005) 191 [astroph/0412109].
 [42] M. V. Kossov, Approximation of photonuclear interaction crosssections, European Physical Journal A 14 (2002) 377.
 [43] M. B. Chadwick et al., Handbook on photonuclear data for applications crosssections and spectra, IAEATECDOC117, wwwnds.iaea.org/photonuclear/
 [44] R. Alves Batista, D. Boncioli, A. di Matteo, A. van Vliet and D. Walz, Effects of uncertainties in simulations of extragalactic UHECR propagation, using CRPropa and SimProp, JCAP 1510 (2015) 10 063 [arXiv:1508.01824].
 [45] A. M. Hopkins and J. F. Beacom, On the normalisation of the cosmic star formation history Astrophys. J. 651 (2006) 142 [arXiv:astroph/0601463]
 [46] L. A. Anchordoqui, J. F. Beacom, H. Goldberg, S. PalomaresRuiz and T. J. Weiler, TeV gammarays from photodisintegration/deexcitation of cosmicray nuclei, Phys. Rev. Lett. 98 (2007) 121101 [astroph/0611580].
 [47] M. Settimo and M. De Domenico, Propagation of extragalactic photons at ultrahigh energy with the code, Astropart. Phys. 62 (2015) 92 [arXiv:1311.6140].
 [48] K. Dolag, D. Grasso, V. Springel and I. Tkachev, Constrained simulations of the magnetic field in the local Universe and the propagation of UHECRs, JCAP 0501 (2005) 009 [astroph/0410419].
 [49] G. Sigl, F. Miniati, and T. A. Ensslin, Ultrahigh energy cosmic ray probes of large scale structure and magnetic fields, Phys. Rev. D70 (2004) 043007 [astroph/0401084].
 [50] M. Lemoine, Extragalactic magnetic fields and the second knee in the cosmicray spectrum Phys. Rev. D 71 (2005) 083007 [astroph/0411173].
 [51] K. Kotera and M. Lemoine, Phys. Rev. D 77 (2008) 023005 [arXiv:0706.1891].
 [52] S. Mollerach and E. Roulet, Magnetic diffusion effects on the ultrahigh energy cosmic ray spectrum and composition, JCAP 1310 (2013) 013 [arXiv:1305.6519].
 [53] R. Alves Batista and G. Sigl, Diffusion of cosmic rays at EeV energies in inhomogeneous extragalactic magnetic fields, JCAP 1411 (2014) 031 [arXiv:1407.6150].
 [54] A. J. Koning and S. Goriely, private communication.
 [55] T. Belgya et al., Handbook for calculations of nuclear reaction data, RIPL2. IAEATECDOC1506, wwwnds.iaea.org/RIPL2