CRPropa 3 - a Public Astrophysical Simulation Framework for Propagating Extraterrestrial Ultra-High Energy Particles

We present the simulation framework CRPropa version 3 designed for efficient development of astrophysical predictions for ultra-high 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 high-performance 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.


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 cut-off above ∼ 40 EeV (1 EeV = 10 18 eV) [1]. Cosmic ray arrival directions appear to be rather isotropic with a few exceptions only. As examples we mention a dipole signal above 8 EeV [2], and a hot spot observed in the Northern hemisphere for energies E > 57 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 ∼ 4 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 (1 PeV = 10 15 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 ∼ 1 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].

JCAP05(2016)038
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 ∼ 5 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 one-dimensional 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 "cut-off" of the spectrum may not be dominated by interactions with background fields as predicted by the so-called 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 [10,15,16]. 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 magneto-hydrodynamics [17,18]. Injection spectra at sources are not yet established, but can be obtained, e.g., from shock acceleration theory.
Furthermore, at energies around 1 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 three-dimensional 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 "four-dimensional 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 ultra-high energy cosmic ray propagation in CRPropa 3 are demonstrated in section 4 in a few example applications. Finally, results are summarized in section 5.

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, electron-positron pairs and neutrinos. Included interactions are pair production, photo-pion 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 one-dimensional (1D) mode and in three-dimensional (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 e + e − 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 e + e − 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

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 user-defined 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 user-provided 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 multi-core processors can therefore be adequately utilized by just running multiple simulation instances in parallel. However, CRPropa 3 enables shared-memory multiprocessing using OpenMP 1 for easy multi-core 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]).

JCAP05(2016)038
General magnetic fields Uniform Magnetic field is a position independent single vector magnetic field.

Grid
Provides a periodic magnetic field grid with trilinear interpolation, equal spacing, and different sizes along each axis.

Modulated Grid
Modulates a large scale vector field by a periodic small scale scalar field. Turbulent A random magnetic field with a turbulent spectrum [23].

Logarithmic Spiral
Magnetic field model of axisymmetric (ASS) or bisymmetric (BSS) logarithmic spiral shape.

Pshirkov 2011
Pshirkov et al. magnetic field model, consisting of a large-scale regular (disk and halo) field [16]. The axisymmetric (ASS) and the bisymmetric (BSS) disk model can be chosen.

JF 2012
Implementation of the Jansson & Farrar magnetic field model, consisting of a large-scale regular and random (striated) field and a smallscale random (turbulent) field [10,15]. Table 1. Implemented general and galactic magnetic field types.
CRPropa 3 is written in C++ and interfaced to Python using SWIG. 2 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 cross-language 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.

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 uses R-Trees for fast access to RAMSES [26,27] magnetohydrodynamical simulations data. Quimby 4 [28] is a multi-resolution grid library for fast access to huge compressed JCAP05(2016)038 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 B com (z) = B(0)(1 + z) m , where B com is the comoving magnetic field and m 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.

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 E i /Z i are stored in matrices L i which are used to map discrete directions (pixel) indexed with n outside the galaxy to discrete observed directions indexed with m on Earth. The set of matrices {L 1 · · · L N } form the 'galactic lens' that completely describes the deflection of extragalactic cosmic rays in the galactic magnetic field model. Using the matrices, a vector p i eg of the probabilities to observe a cosmic ray at energy E i from direction n at the edge of the galaxy can be transformed by a matrix vector multiplication to obtain the probability distribution p i obs for energy E i . To avoid artificial distortions of the energy spectrum in this approach, all matrices have to be normalized by the maximum of unity norms L i 1 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 tem-5 http://eigen.tuxfamily.org. plate library for linear algebra to minimize memory consumption and disk storage space while maximizing the computational performance for the matrices-vector multiplications.

Photonuclear interactions
Photonuclear interactions with the CMB and IRB are the dominant energy loss processes for ultra-high 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 [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 are 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 = 0.2−200 MeV in logarithmic steps of ∆ log 10 ( /MeV) = 0.01. This is done for all 169 isotopes in the range A = 12 − 56, Z ≤ 26 with a lifetime of τ > 2 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, d, t, 3 He and 4 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 A < 12, the same set of cross sections as in CRPropa 2 is used, with the exception of 6 Li, which is now modeled using the parametrization from Kossov [42]; as well as omitting 6 He and 9 Li, which, due to their short lifetime τ < 2s, exhibit negligible photodisintegration. As a result, TALYS is not used for any isotope with A < 12, 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 photo-hadronic 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 -7 -JCAP05(2016)038 differences in typical UHECR scenarios is at level of 10% in the spectral flux; see ref. [44] for more details.

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 (1 + z) −1 , where z is the redshift whose evolution is given by in the standard ΛCDM cosmology. Here, the Hubble parameter at present time is H 0 ≈ 67.3 km/s/Mpc, Ω m ≈ 0.315 is the density of matter in the universe, encompassing both baryonic and dark matter, and Ω Λ ≈ 0.685 is the cosmological constant, assuming a flat universe (Ω tot = Ω m + Ω Λ = 1) [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 n( , z) = (1 + z) 2 n( /(1 + z), z = 0). 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 non-trivial. Similarly to CRPropa 2, we consider the evolution of the IRB through a global scaling factor s(z), 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

IRB (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 section 4.2.
In the more extreme case of high redshift sources (z ≥ 2) 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 one-dimensional environment this is trivial, since the redshift corresponds to the distance of the source to the observer. In a three-dimensional 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 CR-Propa 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 -8 -

JCAP05(2016)038
extends the notion of a three-dimensional observer to four dimensions, so that the observer is considered a hypervolume composed by a sphere of a given radius R obs and a redshift window of size ∆z obs around z = 0. The calculations for z < 0 is obtained by extrapolating the corresponding quantities down to the desired value.
The redshift window (∆z obs , corresponding to ∆t obs ) 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 t H , then ∆t obs /t H ∼ 0.1 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 (E 10 EeV) we are limited by interaction horizons, which are smaller than 100 Mpc, corresponding to a small redshift (z 0.02). For E 1 EeV 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 ∆z obs in a scenario with null magnetic field. The differences between ∆z obs = 0.05 and ∆z obs = 0.20 are of the order of 10 −5 at 1 EeV and 10 −6 at 0.1 EeV.
CRPropa 3 uses comoving coordinates internally. This implies that the flux dilution of comoving magnetic fields with redshift ∝ (1 + z) 2 is implicitly taken into account. Additionally, a redshift scaling of the form B com (z) = B(0)(1 + z) m can be applied to any magnetic field to account for a general field damping resulting in a total magnetic field evolution of B(z) = B com (z)(1 + z) 2 . 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 f (z) ∝ (1 + z) m , where e.g. the star formation rate for z < 1 is typically described with m = 3.4 [45]. Note, that the evolution parameter m for the source density and magnetic field need not be the same, and are completely unrelated.

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 photo-pion 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 electron-positron 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 γ b via pair production γγ b → e + e − and double pair production γγ b → e + e − e + e − . The electrons lose energy via inverse Compton scattering e ± γ b → e ± γ, triplet pair production e ± γ b → e ± e + e − and synchrotron radiation from gyrating in magnetic fields.

JCAP05(2016)038
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 3.6 compared to the baseline in a benchmark application. The speed up was achieved by replacing the on-the-fly 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 onthe-fly calculation and the interpolation is smaller than 10 −8 and thus negligible.
To access the cascade codes, the secondary photons and electrons generated in ultra-high 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.

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 γ = −1 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 photo-meson 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 DINT-EleCa 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 Iron ν from β decay ν from π ± decay All ν UHECR Figure 4. Histograms of the observed neutrino flux from proton (left) and iron (right) primaries separated by individual creation processes. As comparison, the propagated UHECR spectrum is show as a green histogram.
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.

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 large-scale matter density. In this scenario the sources are distributed randomly with a source density of 10 −3 Mpc −3 following the large-scale structure formation simulations from Dolag et al. [48], which have been constrained by the observed large-scale density field with a radius of 110 Mpc around the Milky Way.

JCAP05(2016)038
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 Coma-like galaxy cluster is of the same order as indicated by Faraday rotation measures. To emphasize cosmic-ray 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 large-scale 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 256 3 cells covering a volume of ∼ (132 Mpc) 3 . A higher-resolution magnetic field with Fourier modes taken from a Gaussian distribution with |B(k)| 2 given by a Kolomogov power spectrum and random polarization, with a coherence length of 500 kpc and a total volume of ∼ (13.2 Mpc) 3 , 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 time-independent 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: with N the number of injected particles and E the particle energy. For this example the particles are injected with a spectral index of γ = −1.5 and cutoff energy E cut = 780 EeV down to a minimum energy of E 0 = 1 EeV. We consider the observer to be a sphere with a given radius R obs , 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 R obs 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 trade-off is to take R obs s , where s 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 s ∼ 300 kpc. For this reason we set R obs = 100 kpc.
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 large-scale and turbulent small-scale 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 coordinate-independent angular power spectra, with multipoles a lm for these sky maps are shown in figure 6, normalized to C 0 . The sky maps and power spectra given here are for the full energy range (≥ 1 EeV) and, as the flux is decreasing with energy, are dominated by cosmic rays with energies 10 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 r 1 = 3 (C 1 /C 0 ) (see the first bin of figure 6 for C 1 /C 0 ) for the proton injection scenario before (after) deflections in the GMF is r p 1 = 0.062 (0.066) and for the iron injection scenario is r Fe 1 = 0.067 (0.069). For comparison the expected 99% confidence level upper bounds that would result from fluctuations of isotropic distributions with 10 3 , 10 4 and 10 5 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).

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 ∆z obs = 0.1, i.e. only events that are recorded in the redshift range −0.1 ≤ z ≤ 0.1 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 z = 2. The spectra for both scenarios are shown in figure 7, and the mean and variance of the logarithm of the mass composition, ln A and σ(ln A), observed at Earth for the iron cases are shown in figure 8, together with the universal spectrum, i.e., the spectrum obtained from a one-dimensional 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 -16 -JCAP05(2016)038 E ≈ 10 19.1 eV is due to kinematics: an iron nucleus injected with an energy E inj of several hundred EeV typically suffers complete photodisintegration into 56 nucleons of energy E inj /56 each. Due to the hard source spectrum ∝ E −1.5 these secondary protons dominate the observed flux up to the energy 780 EeV/56 ≈ 10 19.1 eV that corresponds to the cutoff in the source spectrum of the initial iron nuclei. One can see that for E 10 19 eV in the 3D case ln A 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 σ(ln A).
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 E 10 19 eV, the 1D scenario has a slightly larger ln A 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 (E ∼ 1 EeV) 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 B com (z) = B(0)(1+z) m , as discussed in section 3.5. We have not considered the redshift evolution of the field, which is equivalent to taking m = 0 in the aforementioned equation. The parameter m is expected to be non-zero if the magnetic field is significantly affected by magnetohydrodynamic processes taking place during structure formation. Moreover, this approximation is adequate for very high energies (E 10 19 eV), 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.

Summary
The simulation of galactic and extragalactic cosmic ray propagation plays an essential role in understanding astrophysical processes at ultra-high 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 large-scale 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 cosmo-

A Giant dipole resonance parameters
Photodisintegration is the most important interaction for cosmic ray nuclei with energies E > 10 19 eV. Cross sections for this interaction are dominated by the giant dipole resonance for photons with energies < 30 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 RIPL-2 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 12 C and 28 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. 28 Si TALYS 1.6 (adjusted) TALYS 1.6 (default) TALYS 1.0 (default) Kossov IAEA atlas Figure 9. Comparison of total photodisintegration cross sections for 12 C (left) and 28 Si (right) with the evaluated experimental data compiled in the IAEA atlas [43]. TALYS 1.0 (default) and TALYS 1.6 (adjusted) correspond to models implemented in CRPropa 2 and CRPropa 3, respectively. Alternatively, the Kossov parametrization [42] can be used.  Table 2. Giant dipole resonance parameters used with TALYS (as parameters for the Kopecky-Uhl generalized Lorentzian model of the E1-strength function): peak energy E i , peak cross section σ i and width Γ i for resonances with a single (i = 0) or a split peak (i = 0, 1). Default values from the RIPL-2 database [55] are replaced, if available, with the total cross section parameters from the atlas of giant dipole resonance parameters. Note that for isotopes not listed, as well as for higher order contributions, TALYS uses a compilation of formulas listed in [40].