text:
\section{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 $\sim 40$~EeV ($1~{\rm EeV}= 10^{18}\,$eV) \cite{Abraham:2010mj}. Cosmic ray arrival directions appear to be rather isotropic with a few exceptions only. As examples we mention a dipole signal above $8$~EeV~\cite{ThePierreAuger:2014nja}, and a hot spot observed in the Northern hemisphere for energies $E>57$~EeV~\cite{Abbasi:2014lda}. 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 $\sim 4$~EeV~\cite{Aab:2014aea}.

In recent years, high energy neutrinos have also been observed with a flux well above expectations from atmospheric showers \cite{Aartsen:2014gkd}. Such extraterrestrial neutrinos with energies in the PeV regime ($1~{\rm 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 \cite{Aartsen:2014gkd}, or with cosmic ray arrival directions \cite{Aartsen:2015dml} 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 $\sim 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~\cite{Giacinti:2011ww,Abreu:2012lva}. Taking advantage of the isotropic arrival distributions, at least bounds on the density of sources were derived depending on the cosmic ray energy~\cite{Abreu:2013kif}.

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 $\sim5$~EeV \cite{Jansson2012}, 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~\cite{Greisen1966,Zatsepin1966} but rather be caused by the limited source energies \cite{Aab:2015bza}.

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.\ \cite{Gilmore:2011ks}. 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 \cite{Pshirkov11,Jansson2012,Jansson2012a}. 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 \cite{Dolag:2008ya,Miniati:2011}. 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 \cite{Agashe:2014kda}.

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~\cite{Kampert2012} 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 \ref{Sec:InhFeat} we will briefly recap the inherited functionality and physics from the previous version CRPropa 2. Section \ref{Sec:NewFeat} 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 \ref{Sec:Examples} in a few example applications. Finally, results are summarized in section \ref{Sec:Summary}.

\section{Inherited Features from CRPropa 2}
\label{Sec:InhFeat}
In the following we will introduce the features of CRPropa version 2~\cite{Kampert2012} 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~\cite{Lee1998}. All these features of CRPropa 2 are available in CRPropa 3 as well.

\section{New Features in CRPropa 3}
\label{Sec:NewFeat}
\subsection{Code structure and steering}
\label{Sec:CodeStruc}
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 time-dependencies as input.

A graphical visualization of the propagation process is given in figure~\ref{fig:visu}. 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.

\begin{figure}[tb]
\begin{centering}
\includegraphics[scale=0.5]{fig1.pdf}
\caption{Illustration of the CRPropa 3 modular structure. Each module contained in the module list acts on the candidate class. The \texttt{isActive} flag serves as break condition and is checked after each cycle of the module list.}
\label{fig:visu}
\end{centering}
\end{figure}

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\footnote{\url{http://openmp.org}} 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~\cite{Mucke2000}).

CRPropa 3 is written in C++ and interfaced to Python using SWIG\footnote{\url{http://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 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 \texttt{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 \ref{Sec:Examples}.

\subsection{Magnetic fields}\label{sec:BFields}
As part of the restructuring, CRPropa 3 now supports any kind of magnetic field.
The only requirement is the implementation of a \texttt{getField} function in the C++ \texttt{MagneticField} interface.
This allows analytical fields, grid like fields, or fields from complex structures of any scale to be included. Table \ref{table:MagFields} summarizes the generic general and galactic magnetic field types implemented in CRPropa 3.
\begin{table}[t!]
\centering
\makebox[0pt][c]{\parbox{1.05\textwidth}{%
    \begin{minipage}[c][4.2cm][t]{0.485\hsize}\centering
    \begin{scriptsize}
        \begin{tabular}{| p{1.3cm}p{5.8cm}  | }
           \hline
           \multicolumn{2}{|c|}{\textbf{General magnetic fields}}
           \\ \hline
           \textit{Uniform}    & Magnetic field is a position independent single vector magnetic field.\\
           \textit{Grid}           & Provides a periodic magnetic field grid with trilinear interpolation, equal spacing, and different sizes along each axis.\\
           \textit{Modulated Grid} & Modulates a large scale vector field by a periodic small scale scalar field.   \\
           \textit{Turbulent}      & A random magnetic field with a turbulent spectrum~\cite{Giacinti11}.   \\
           & \\
           & \\
           & \\ \hline
        \end{tabular}
        \end{scriptsize}
    \end{minipage}
    \hfill
    \begin{minipage}[c][4.2cm][t]{0.52\hsize}\centering
    \begin{scriptsize}
        \begin{tabular}{| p{1.3cm}p{6.115cm}  | }
           \hline
           \multicolumn{2}{|c|}{\textbf{Galactic magnetic fields}}
           \\ \hline
            \textit{Toroidal Halo}     & Toroidal halo field model adopted from~\cite{Prouza:2003yf, Sun:2007mx}.\\
           \textit{Logarithmic Spiral} & Magnetic field model of axisymmetric (ASS) or bisymmetric (BSS) logarithmic spiral shape. \\
            \textit{Pshirkov 2011}          &  Pshirkov et al.\ magnetic field model, consisting of a large-scale regular (disk and halo) field~\cite{Pshirkov11}. The axisymmetric (ASS) and the bisymmetric (BSS) disk model can be chosen. \\
           \textit{JF 2012}   & Implementation of the Jansson \& Farrar magnetic field model, consisting of a large-scale regular and random (striated) field and a small-scale random (turbulent) field~\cite{Jansson2012, Jansson2012a}.   \\  \hline
        \end{tabular}
        \end{scriptsize}
    \end{minipage}
    \hfill
        \caption{Implemented general and galactic magnetic field types.}
    \label{table:MagFields}
}}
\end{table}

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\footnote{\url{https://github.com/rafaelab/saga}} uses R-Trees for fast access to RAMSES \cite{Teyssier:2006us,Fromang:2006aw} magnetohydrodynamical simulations data.
Quimby\footnote{\url{https://forge.physik.rwth-aachen.de/projects/quimby}}~\cite{Muller:2015gkr} is a multi-resolution grid library for fast access to huge compressed magnetic field grids and provides access to smooth particles from the GADGET~\cite{Springel:2000yr,Springel:2005mi} 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_{\rm com}(z) = B(0) (1+z)^{m}$, where $B_{\rm com}$ is the comoving magnetic field and $m$ the magnetohydrodynamic amplification of the field, cf.\ section \ref{Sec.CosEffects}.
Multiple magnetic fields can also be grouped together and their magnetic fields can be superpositioned (added) in a \textit{list}.

\subsection{Galactic propagation}\label{sec:GalProp}
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~\cite{Jansson2012, Jansson2012a} and Pshirkov et al.~\cite{Pshirkov11} 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~\cite{Bretz2014} 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 $\matrix{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 $\lbrace
\matrix{L}_1 \cdots \matrix{L}_N \rbrace$ form the `galactic lens' that
completely describes the deflection of extragalactic cosmic rays in the galactic
magnetic field model. Using the matrices, a vector $\vecN{p}_{\rm eg}^i$ 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
\begin{equation}
	 \matrix{L}_i \cdot \vecN{p}_{\rm eg}^i = \vecN{p}_{\rm obs}^i
	\label{eq:GalacticLensesEquation}
\end{equation}
to obtain  the probability distribution $\vecN{p}_{\rm obs}^i$ 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
$\Vert \matrix{L}_i \Vert_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~\cite{Gorski2005}. The matrices
are accessed in sparse compressed column major format using the
Eigen\footnote{\url{http://eigen.tuxfamily.org}} template library for
linear algebra to minimize memory consumption and disk storage space while
maximizing the computational performance for the matrices-vector
multiplications.

\subsection{Photonuclear interactions}
\label{sec: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~\cite{Kneiske2004} model that was considered in CRPropa 2, the following IRB models are now available as well: Stecker~2005~\cite{Stecker2005}, Franceschini~2008~\cite{Franceschini2008}, Finke~2010~\cite{Finke10},  Dominguez~2011~\cite{Dominguez11} and Gilmore~2012~\cite{Gilmore:2011ks}.

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~\cite{Kalashev2014}.
Additionally, the following improvements regarding photodisintegration were implemented.
In CRPropa 2 photodisintegration cross sections were mainly obtained with TALYS 1.0~\cite{Goriely2008}; see ref.~\cite{Kampert2012} for a detailed description.
These cross sections are updated using the more recent TALYS version 1.6 (cf.\ ref.~\cite{TALYSmanual} 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 \cite{Khan:2004nd}.
The complete list of used parameters is given in appendix~\ref{appendix:GDR}, table \ref{tab:GDRparameters}.
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 $\epsilon' = 0.2 - 200$\,MeV in logarithmic steps of $\Delta \log_{10}(\epsilon'/\mathrm{MeV}) = 0.01$.
This is done for all 169 isotopes in the range $A = 12 - 56$, $Z \leq 26$ with a lifetime of $\tau > 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~\cite{Kossov2002}; as well as omitting $^6$He and $^9$Li, which, due to their short lifetime $\tau < 2s$, exhibit negligible photodisintegration.
As a result, TALYS is not used for any isotope with $A < 12$, as is recommended in ref.~\cite{TALYSmanual}.
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~\cite{Kossov2002}.
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 \cite{iaea-tecdoc-1178} 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.~\cite{SimProp_CRPropa} for more details.

\subsection{Cosmological effects}
\label{Sec.CosEffects}
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
\begin{equation}
 \left| \frac{{\rm d}t}{{\rm d}z} \right| = \frac{1}{H_0(1+z) \sqrt{\Omega_m(1+z)^3+\Omega_\Lambda}}
\end{equation}
in the standard $\Lambda$CDM cosmology. Here, the Hubble parameter at present time is $H_0\approx67.3$ $\mathrm{km/s/Mpc}$, $\Omega_m\approx0.315$ is the density of matter in the universe, encompassing both baryonic and dark matter, and $\Omega_\Lambda\approx0.685$ is the cosmological constant, assuming a flat universe ($\Omega_{\rm tot} = \Omega_{m} + \Omega_{\Lambda} =1  $)~\cite{Agashe:2014kda}.

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(\epsilon,z) = (1+z)^2 \, n(\epsilon/(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 \ref{sec:Interactions}) 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
\begin{align}
  \lambda^{-1}(E,z) &= s(z) (1+z)^3 \lambda^{-1}(E(1+z), z=0) \\
  s(z) &= \begin{cases} \quad\quad 1 & \text{CMB} \\ \frac{\int n(\epsilon,z) {\rm d}\epsilon}{\int n(\epsilon,0) {\rm d}\epsilon} & \text{ IRB} \end{cases}
\end{align}
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. \ref{sec:Benchmark_example}. In the more extreme case of high redshift sources ($z\geq 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 {\it 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 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 three-dimensional observer to four dimensions, so that the observer is considered a hypervolume composed by a sphere of a given radius $R_{\rm obs}$ and a redshift window of size $\Delta z_{\rm obs}$ around $z=0$. The calculations for $z<0$ is obtained by extrapolating the corresponding quantities down to the desired value.

The redshift window ($\Delta z_{\rm obs}$, corresponding to $\Delta t_{\rm 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 $\Delta t_{\rm obs} / t_H \sim 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 \gtrsim 10 \text{ EeV}$) we are limited by interaction horizons, which are smaller than $\lesssim 100 \text{ Mpc}$, corresponding to a small redshift ($z \lesssim 0.02$). For $E\lesssim 1 \text{ 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 $\Delta z_{\rm obs}$ in a scenario with null magnetic field. The differences between $\Delta z_{\rm obs}=0.05$ and $\Delta z_{\rm 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 $\propto (1+z)^2$ is implicitly taken into account. Additionally, a redshift scaling of the form $B_{\rm 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_{\rm 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 \ref{sec:BFields}.

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) \propto (1+z)^{m}$, where e.g. the star formation rate for $z < 1$ is typically described with $m=3.4$~\cite{Hopkins:2006bw}. Note, that the evolution parameter $m$ for the source density and magnetic field need not be the same, and are completely unrelated.

\subsection{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~\cite{Anchordoqui:2006pd}. 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 $\gamma_b$ via pair
production $\gamma\gamma_b \rightarrow e^+e^-$ and double pair production
$\gamma\gamma_b \rightarrow e^+e^-e^+e^-$. The electrons lose energy via
inverse Compton scattering $e^\pm\gamma_b \rightarrow e^\pm\gamma$, triplet pair
production $e^\pm\gamma_b\rightarrow e^\pm e^+e^-$ 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~\cite{Lee1998}, which was also
part of the previous versions of CRPropa, and as a new feature also to
EleCa~\cite{Settimo2013}. DINT calculates the observed spectra based on the transport equations given in ref.~\cite{Lee1998}. 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 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
on-the-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 \ref{Sec.AppSecondary}.

 \section{Example Applications Focussing on New Features}
\label{Sec:Examples}
\subsection{1D simulation including secondary particles}
\label{Sec.AppSecondary}
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 $\gamma = -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 \ref{fig:injected_secondaries} separated into
generating processes and type of primary.
\begin{figure}[tb]
	\includegraphics[width=.49\textwidth]{fig2.pdf}
	\includegraphics[width=.49\textwidth]{fig3.pdf}
	\caption{Histograms of the number of secondaries injected into the
	electromagnetic cascade from the propagation of proton (left) and iron
(right) primaries for the individual creation processes integrated over the
whole propagation distance, see section \ref{Sec.AppSecondary} for details. As comparison, the injected primary spectrum is shown as well (green). The observed photon and neutrino fluxes from this simulation are shown in figures \ref{fig:observed_photons} and \ref{fig:observed_neutrino_flux}.}
	\label{fig:injected_secondaries}
\end{figure}
For protons, the majority of secondary particles are low energy electrons
created in electron pair production. Higher energy electrons are also generated
in the $\beta$-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 $\beta$-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~\ref{fig:observed_photons} as histograms for the combined
propagation and gray surface for the DINT only propagation.
\begin{figure}[tb]
	\includegraphics[width=.49\textwidth]{fig4.pdf}
	\includegraphics[width=.49\textwidth]{fig5.pdf}
	\caption{Histograms of the observed photon flux from proton (left)
	and iron (right) primaries separated by individual processes for the combined
propagation. As comparison, the DINT only propagation as in previous CRPropa
versions is shown as a gray background, and the propagated UHECR spectrum is shown as a green histogram.}
	\label{fig:observed_photons}
\end{figure}
For energies above a few hundred TeV ($1~{\rm TeV}= 10^{12}\,$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~\ref{fig:observed_neutrino_flux}. In the case of proton primaries,
approximately four times more neutrinos are generated in decays of $\pi^\pm$
compared to $\beta$-decays. In the case of iron primaries, approximately ten times
more neutrinos are generated from $\beta$ decay than from the decay of
$\pi^\pm$.  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.
\begin{figure}[tb]
	\includegraphics[width=.49\textwidth]{fig6.pdf}
	\includegraphics[width=.49\textwidth]{fig7.pdf}
	\caption{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. }
	\label{fig:observed_neutrino_flux}
\end{figure}

\subsection{3D simulation including extragalactic and galactic deflections}
\label{sec:Benchmark_example}%

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}~\text{Mpc}^{-3}$ following the large-scale structure formation simulations from Dolag et al.~\cite{Dolag2005}, which have been constrained by the observed large-scale 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.~\cite{Miniati2004}. 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~\ref{sec:BFields}) of $256^3$ cells covering a volume of $\sim (132\text{ Mpc})^3$. A higher-resolution magnetic field with Fourier modes taken from a Gaussian distribution with $\langle |B({\bf k})|^2\rangle$ given by a Kolomogov power spectrum and random polarization, with a coherence length of 500~kpc and a total volume of $\sim (13.2\text{ 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~\ref{sec:Interactions}) 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~\cite{Gilmore:2011ks}. 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:
\begin{equation}\label{eq:PowerLawInjection}
	\frac{\text{d}N}{\text{d}E} \propto \left\{
	\begin{array}{l l}
 	(E/E_0)^{\gamma} & E \leq E_{\text{cut}} \\
 	(E/E_0)^{\gamma}\exp(1-E/E_{\text{cut}}) & E > E_{\text{cut}}
	\end{array}
	\right.
\end{equation}
with $N$ the number of injected particles and $E$ the particle energy. For this example the particles are injected with a spectral index of $\gamma = -1.5$ and cutoff energy $E_{\text{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_{\rm 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_{\rm 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_{\rm obs} \lesssim \ell_s$, where $\ell_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 $\ell_s \sim 300 \text{ kpc}$. For this reason we set $R_{\rm obs} = 100 \text{ 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~\ref{sec:GalProp} has been used. The implemented galactic magnetic field model is the Jansson and Farrar 2012 model~\cite{Jansson2012, Jansson2012a} 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 \ref{fig:3D_skymaps}. 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.

\begin{figure}[tb]
	\subfigure[3D proton, no galactic magnetic field]{
    	\includegraphics[width=.49\textwidth]{fig8.pdf}
  	}
  	\subfigure[3D proton, with galactic magnetic field]{
    	\includegraphics[width=.49\textwidth]{fig9.pdf}
  	}
  	\subfigure[3D iron, no galactic magnetic field]{
    	\includegraphics[width=.49\textwidth]{fig10.pdf}
  	}
  	\subfigure[3D iron, with galactic magnetic field]{
    	\includegraphics[width=.49\textwidth]{fig11.pdf}
  	}
	\caption{Sky maps in galactic coordinates for pure proton injection (top) and pure iron injection (bottom) before (left) and after (right) deflections in the galactic magnetic field are taken into account. All sky maps have been normalized so that the bin with the maximum number of hits is set to one. For better visibility we restrict the color scale to values between 0.3 and 0.7. The simulations use a spectral index of $\gamma = -1.5$ and cutoff energy $E_{\text{cut}} = 780$~EeV down to a minimum energy of $E_0 = 1$~EeV. See text for further details.}
	\label{fig:3D_skymaps}
\end{figure}

The coordinate-independent angular power spectra,
\begin{equation}
	C_l = \frac{1}{2l+1}\sum^l_{m=-l}|a_{lm}|^2
\end{equation}
with multipoles $a_{lm}$ for these sky maps are shown in figure\ \ref{fig:3D_powerspectra}, normalized to $C_0$. The sky maps and power spectra given here are for the full energy range ($\geq1$~EeV) and, as the flux is decreasing with energy, are dominated by cosmic rays with energies $\lesssim 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\sqrt(C_1/C_0)$ (see the first bin of figure\ \ref{fig:3D_powerspectra} for $C_1/C_0$) for the proton injection scenario before (after) deflections in the GMF is $r_1^{\text{p}} = 0.062$ (0.066) and for the iron injection scenario is $r_1^{\text{Fe}} = 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 \ref{fig:4Dspec} and \ref{fig:4Dcomp} together with the spectra and compositions for the 4D cases (cf.\ section \ref{sec:4d}).

\begin{figure}[tb]
	\centering
	\includegraphics[width=.6\textwidth]{fig12.pdf}
	\caption{Angular power spectra of the cosmic ray arrival directions for pure proton (green) and pure iron (orange) injection at the sources without (dashed lines) and with (solid lines) deflections in the galactic magnetic field taken into account. Additionally, the $99\%$ confidence level upper bounds that would result from fluctuations of isotropic distributions with $10^3$, $10^4$ and $10^5$ events are shown (black solid lines). The simulations are done using a spectral index of $\gamma = -1.5$ and cutoff energy $E_{\text{cut}} = 780$~EeV down to a minimum energy of $E_0 = 1$~EeV. See text for further details.}
	\label{fig:3D_powerspectra}
\end{figure}

\subsection{Cosmological effects using 4D simulations}
\label{sec:4d}
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 $\Delta z_{\rm obs} = 0.1$, i.e.\ only events that are recorded in the redshift range $-0.1 \leq z \leq 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~\ref{fig:4Dspec}, and the mean and variance of the logarithm of the mass composition, $\langle \ln A \rangle$ and $\sigma(\ln A)$, observed at Earth for the iron cases are shown in figure~\ref{fig:4Dcomp}, together with the universal spectrum, i.e., the spectrum obtained from a one-dimensional simulation assuming a uniform distribution of sources.
\begin{figure}[tb]
	\centering
	\includegraphics[width=0.80\columnwidth]{fig13.pdf}
	\caption{Simulated spectra for the proton (lower panel) and iron (upper panel) injection. The orange solid lines correspond to the 4D simulation, the 3D case is represented by green dashed lines, and dot-dashed black lines correspond to the 1D spectrum. The spectra are normalized at 10$^{20} \text{ eV}$. Parameters used in these simulation are: $\gamma=-1.5$, $E_{\rm cut}=780\text{ EeV}$, and $E_0 = 1$~EeV.}
	\label{fig:4Dspec}
\end{figure}
\begin{figure}[tb]
	\centering
	\includegraphics[width=0.49\columnwidth]{fig14.pdf}
	\includegraphics[width=0.49\columnwidth]{fig15.pdf}
	\caption{Mass composition observed at Earth for the case of a source emitting only iron. The average of the logarithm of the mass number, $\langle \ln A\rangle$, and its standard deviation, $\sigma(\ln A)$, are shown in the left and right panels, respectively. Orange solid lines correspond to the 4D simulation, green dashed lines correspond to the 3D case, and dot-dashed black lines represent the 1D spectrum. Parameters used in these simulation are: $\gamma=-1.5$, $E_{\rm cut}=780\text{ EeV}$, and $E_0 = 1$~EeV.}
	\label{fig:4Dcomp}
\end{figure}

In figure~\ref{fig:4Dspec} 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~\ref{fig:4Dcomp}. Note that the common feature of an abrupt increase in average mass at around $E \approx 10^{19.1}\,\text{eV}$ is due to kinematics: An iron nucleus injected with an energy $E_{\rm inj}$ of several hundred EeV typically suffers complete photodisintegration into 56 nucleons of energy $E_{\rm inj}/56$ each.
Due to the hard source spectrum $\propto E^{-1.5}$ these secondary protons dominate the observed flux up to the energy $780\,\text{EeV}/56 \approx 10^{19.1}\,\text{eV}$ that corresponds to the cutoff in the source spectrum of the initial iron nuclei.
One can see that for $E \lesssim 10^{19}\text{ eV}$ in the 3D case $\langle \ln A \rangle$ 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 $\sigma(\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 \lesssim 10^{19}\text{ eV}$, the 1D scenario has a slightly larger $\langle \ln A \rangle$ 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 \sim 1 \text{ 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~\cite{Lemoine:2004uw,Kotera:2007ca,Mollerach:2013dza,Batista:2014xza}.

Although we have taken redshift effects into account, the comoving magnetic field can also evolve with redshift approximately as $B_{\rm com}(z) = B(0) (1+z)^{m}$, as discussed in section~\ref{Sec.CosEffects}. 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 \gtrsim 10^{19} \text{ 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.

\section{Summary}
\label{Sec: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 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 \textit{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 multi-messenger 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 $\sim 100$~TeV.

Information on downloading the code, usage and example applications can be found at \url{http://github.com/CRPropa/CRPropa3}. Questions and comments can be submitted to the ticketing system \url{https://github.com/CRPropa/CRPropa3/issues}. The features of CRPropa 3 can also be explored online without installation at \url{https://vispa.physik.rwth-aachen.de}.

redpajama_set_name: RedPajamaArXiv


Over the last several years, Hopsin has been offering fans regular glimpses into his headspace via an ongoing series of video singles entitled Ill Mind of Hopsin. Recently, the Funk Volume rapper/producer's been lost in bittersweet memories of an Old Friend who passed away due to a drug overdose. On installment number six in the aforementioned series, the artist revisits the painful experience of watching his homie self-destruct in slow motion, flowing over his own mournful string loops. Director George Orozco provides a visual accompaniment to the artist's bracing narrative bars. Out now on iTunes, this joint isn't attached to any larger project, but fans should keep their browsers locked for new song releases leading up to Hopsin's Knock Madness LP, which is scheduled for release this November.


Ksenia judges books by their covers – and then she elevates them into works of art.
Hand-crafted, one of a kind, leather-bound novels line the walls of Ksenia Kopystyńska's downtown limestone brick studio in Kingston, Ontario. The towering windows of Ars Libri Studio (latin for "the art of books") feature her handiwork, giving passers-by the opportunity to browse beautiful hand-bound books and other leathercraft items. Even the exterior walls of this heritage building seem to breathe a legacy almost forgotten, underlining an important and historical connection to the work she does daily.
Ksenia began her career almost 30 years ago with a master's degree in Book Conservation and Binding and taught this professionally in Canada for many of these. A sudden but lengthy trip back to Europe to attend to family matters, and a period of unemployment on her return to Canada, meant that she had to essentially restart her career from scratch. Through the EI system, Ksenia discovered Momentum's Self employment Program.
Her works of art are more than simple dust covers. Ksenia prides herself in working with many of the most prolific Canadian and American book collectors and academic institutions. The complexity that goes into every piece she makes means that, unlike the business plans of many entrepreneurs, hers do not include any further expansion.
Ksenia's future goals do, however, include being able to work with the most prominent collections within Canada, and with the collectors who appreciate the value of these significant books and documents.
To launch her business, Ksenia was offered a Micro Business Loan from Momentum. She freely admits that costs of launching a business can be substantial, and that the real value of the program and microloan is in the confidence they give. "It's nice to know that somebody believes in you," she says.
With an initial investment of $30,000 from a friend, and the microloan offered by Momentum for some basic materials, Ksenia moved to Kingston where she determined the best market for her wares would be found. Her market research was accurate – over 360 people attended her business launch Open House.
Ksenia's plans include expanding the brand to offer more hand-stitched designer leather bags and purses. But currently these items are only offered to high end clients seeking specialty creations.
To see more of Ksenia's works of art, check out her business website, Ars Libri Studio.
Like Ksenia, do you have a passion to produce something beautiful?
Get in touch with our Business Development department to discover SELF-EMPLOYMENT today!


Does it make it better that it did the kicking with a fabulous pair of heels? Sure, why not. After making its debut at the Steven Universe SDCC panel, "Legs From Here to Homeworld" hit the Cartoon Network app last night for public consumption. In typical Steven Universe fashion, it managed to pack a ludicrous amount of content into 11 minutes. Let's get into it, because we have a lot to talk about.
The episode kicks off with a overwhelmingly happy Blue Diamond snuggling Steven and making everyone cry her tears. Garnet panics about the possibility of Blue hurting him, but Bismuth assures her she won't. Yellow eventually buzz-kills the whole thing to start asking why Pink faked her shattering and looks like she does. Steven does his best to explain moms and how he's not fully Pink, but Yellow assumes they can just bring back her memories and such later. Yellow Mom is nothing if not stubborn.
Blue wonders how they survived the corruption attack, which it turns out was meant to destroy, not corrupt. Steven brings them to Centipeetles ship and brings her out, hoping they can help reverse the damage done. Yellow and Blue try to heal her at Steven's behest. They manage to cure her form and ability to communicate, but her mind is still stuck in a loop from the evacuation. It's not until Steven adds his powers to the other Diamonds that her mind clears and Nephrite emerges.
Yes, theory confirmed! Centi is a Nephrite, and she's amazing! She apologizes for failing to evacuate but assumes that since the Diamonds are on Earth, they won the war. The moment they all stop touching her, though, Nephrite corrupts again. Steven wants to heal her again, and says healing all the corruption is the Diamonds' responsibility.
Yellow says they can't heal them without all 4 Diamonds. Yes, it happened! The White Diamond name drop! Pearl immediately freaks, saying White is nothing like other gems, or even the other Diamonds. Blue and Yellow aren't very optimistic either. Steven insists they go to Homeworld to talk to her, though. When Yellow questions how, since their arm ships were destroyed, Steven takes them out to the desert to uncover Pink's leg ship.
Look, we all knew it was a leg ship. We all expected it. And oh my god is it magnificent.
They fly the ship back to the beach, where the Fellowship forms. All the Crystal Gems except for Bismuth will go. She decides to stay behind to watch over Lapis and Peridot's gems. Lion stays, too. Connie says goodbye to her parents and Steven says goodbye to Greg. Garnet does as well, handing off Cat Steven for him to take care of. With all the goodbyes finished, it's time to take off! Which requires a gloriously goofy running start. We could not love the leg ship more than we do.
During the trip, Yellow and Blue paint an unflattering picture of White Diamond as a difficult person. They worry how she'll react to the truth of the rebellion on Earth. Yellow even flat out states that only she should be allowed to do the talking. When they arrive at Homeworld, we get our first real look of the planet, which is shattered into three pieces barely held together. And yes, the ship lands as fabulously as it took off. A loud cheering can be heard, and they exit the ship to find gems belonging to every Diamond gathered to meet them.
White Pearl shows up to retrieve "Pink Diamond," and only Pink. For real, White Pearl might be the best gem design yet. She brings Steven to White Diamond, and we get our first look at the Homeworld matriarch. Yes, she's amazing. Absolutely amazing. Her design plain works for us, and Bo was admittedly a bit skeptical just seeing the grainy shots of her from SDCC.
White doesn't act very angry but rather…dismissive. After a quick, chiding exchange welcoming Pink home, Steven's taken away to a chamber with a view of the White Diamond statue/ship/whatever, and we're back on hiatus.
Garnet asked Greg to water Cat Steven for her. So, so adorable.
Speaking of Cat Stevens, check that possible foreshadowing for White Pearl. Or was fake Connie in Rose's room foreshadowing? It's all foreshadowing!
Yellow Diamond was trying and barely failing to move her ship telekinetically. This is big news. The Diamonds can move their ships with their freaking minds.
The Diamonds used the warp pad! We choose to believe they shapeshifted smaller in order to fit in the house and use its warp pad.
Look at all those Peridots and Aquamarines!
We really want to know just how similar Steven acted compared to Pink. The leg ship launch was clearly a "Pink hasn't changed" moment. What else, though? We need details.
Who else thinks White's Pearl was the original Pink Pearl? Between the gem placements and White Pearl's scar, we're thinking our Pearl was originally White's Pearl. The two switched for some reason as yet undetermined. It may even be that our Pearl was created for White but gifted to Pink before she could actually serve White. White/OriginallyPink? Pearl was punished for something, hence the scar.
Alright, for once we're just going to dive right in to the thing everyone is talking about. We saw White Diamond. In the span of 11 minutes we went from not even having a name drop to a Steven conversation with her. A very one-sided conversation, but a conversation nonetheless. Now begins the theorizing.
Namely, what the hell is up with White Diamond and her Pearl?
We can't shake just how…robotic they both were. Everything from the moment Steven arrived on Homeworld felt like a scripted event. It has us wondering deeply about the White Diamond we saw. Her spiel to Steven felt like something she'd said a thousand times to fit numerous little rebellious incidents from Pink. The Gem War was such a momentous event to everyone we've seen, including Blue and Yellow, and yet White acted like "Pink" had just run away for the night. She never even responded to anything Steven said, as if it were a pre-recorded conversation.
Then also consider how immobile she was. Not a single part of her moved besides her face. The close-up of her legs almost resembled a statue. A lot of fans have theorized that this wasn't the real White Diamond at all, but some kind of projection. We suppose that could make sense considering how bright White's head was. We couldn't even see her gem. Everything about her seemed off in some way. Even Blue and Yellow talked about how different she was.
Is this just her personality, as the greatest and strongest of all gems? Has she been affected in some way by some incident, like the corruption attack? Is she hiding for some reason and using this projection?
Obviously there is a White Diamond, and a giant White Diamond at that. We saw her participate in the corruption attack. But there was definitely something off about her, as well as her Pearl. They shared a very similar voice, and White Pearl didn't seem to have much, if any control over herself. Her body never moved the way any of the other gems do. She slid around on the ground like she was on a sledge rather than moving her legs. Like White, only her face actually moved that we saw. Perhaps White controls her? And if White controlled her Pearl, maybe she also controlled the White Diamond we saw. If so, this level of control might tell us what White's powers are. Her Pearl's appearance might also tell us about White's capabilities.
In fact, her eyes may be a hint as to what exactly she sees and how she controls it, as well as her nickname for Pink.
White's eyes resemble a blind person's with the clouded iris and white pupils. Is it possible she does not see gems as we see them, but rather just their aura? If so, that would explain why she does not react at all to Steven's appearance. She may only see whatever aura Pink always had. If she can see auras rather than physical forms, it could also explain a power involving control over a gem's form. It would make sense if the first gem could influence design in order to begin and refine the creation process for those who follow.
Consider White Pearl's form; why does she have a scar? Gems do not have permanent physical forms. They do not scar, not that we've seen. A scar like White Pearl's especially would seem to only happen from the type of injury that would cause a gem to retract their physical form and start reform. That she has not suggests White can control the physical forms of other gems, and possibly by controlling their aura.
This would make sense, because it seems to be the last element missing to fully cure the corrupted gems.
We love the scene curing Nephrite, even if we're sad to see her corrupt again. It showed the combined effort needed to cure these gems, and the way each Diamond controls one aspect of the cure. Yellow was able to restore Nephrite's physical form, which makes sense given her connection to the de-fusing weapons we've seen Jasper and Topaz interact with (both yellow gems). Blue calmed her emotions and ability to express herself, which makes sense given Blue's ability to project her feelings and emotional expression onto other gems. Steven calmed her thoughts, which, again, makes sense given his ability to astral project into other gem's minds and interact with the ways they think about themselves or others. We saw that on display in "Reunited" when Steven encouraged the Crystal Gems via projecting positive thoughts into their minds.
All that was missing in Nephrite's healing was the ability to make it permanent. White being needed to control the permanence of these changes, presumably by fixing a gem's aura, makes a lot of sense.
It also explains what happened in the corruption attack. The combined force of the Diamonds seemed to have destabilized the gems on Earth in each area the three Diamonds had control over: emotional imbalance, form distortion, and a lack of fixed aura. They can't be healed through any one Diamond because they each control a different aspect of what happened to these gems.
What we have developing here is yet another beautiful metaphor for the healing process. When it comes to trauma, there's never just one part of a person affected. Physical trauma changes how people think about themselves and expression their emotions. Emotional trauma impacts how people behave and mental conditioning from abuse can warp people's ability to feel and express it, etc. And all of these kinds of trauma affect someone's 'aura' or 'soul,' that deep part of them that defines who they are. Because damage to the self impacts every area of a person's selfhood and expression—mind, body, spirit, and emotions—healing must occur on all these levels as well.
It's just one more way that Queen Sugar and the Crewniverse interweave therapeutic practices and cognitive behavioral science into a story kids can understand. On a show that embodies relationships as individual characters (the fusions), we're not surprised to see healing embodied as characters. And if we're going to get a Diamond Voltron (and by this point, we're definitely going to see a Diamond Voltron; maybe that's how they heal all the gems? By joining up into a giant robot and blasting Earth with healing powers?), that's one more layer of symbolism: heart, mind, body, and soul all fused (heh) together to make up a the Giantest Woman of them all. Greg might have a heart attack.
It's also interesting to note how the original four Crystal Gems mirror the Diamonds in some way. Garnet as body/form, Amethyst as heart, Pearl as mind, and Rose as soul/aura. Just look at their gem placement—hands, heart, forehead, and belly. Gretchen is intrigued by the seeming difference between Rose's role in the Crystal Gems and what we see Steven doing. It could be that we've misconstrued Steven's role with Nephrite. Maybe he's meant to play the aura role, touching the deep part of someone's soul (which makes sense because Steven) rather than influencing their thoughts directly. Maybe White's role is more mind control than aura control. Perhaps she's controlling White Pearl's mind and only that, which could explain White Pearl's lack of movement (no body control). The more Gretchen things about it, the more she's convincing herself that Steven represents aura.
Or maybe we're both over-thinking this. Then again, it's Steven Universe so probably not. We may not be right, but in terms of analysis, there really is no way to over analyze this show. Everything is foreshadowing after all.
Okay, back to the attack. If White Diamond destabilized auras (or minds) as her part in the attack, then it explains why Yellow and Blue had no idea their attack corrupted anything rather than flat out destroying the gems. So now the question remains about whether White knew. We believe she did. Blue and Yellow suggested she is at the very least known for cruelly lashing out, and also that she feels little to nothing for the lives of the gems under her. If we also consider the demeaning way she spoke to Steven/Pink, it's not exactly a stretch to think that she would have corrupted the auras (again, or minds) of Earth's gems as a "lesson" to Pink. Even if she didn't purposefully do it we doubt she'll care much or want to help heal them.
Of course, we know they will get healed eventually. How? We're not sure quite yet, but we can't wait to see. White's first appearance was incredible and we can't wait to learn more about her either. After all the discussion among the fandom that mainly guessed White would be this cold, unfeeling dictator throwing herself into imperial management to distract from her grief, it's looking like she may be the opposite.
Yellow and Blue talked about White like big sisters warning their little sister about another older, abusive sister. Or perhaps an abusive parent. Their hesitance to admit failure to her, their warnings about how she may react, Yellow's desire to "go first" and take the brunt of White's reaction, it all sounded like fear of immediate, cruel consequences. Even Pink punching a window in "Jungle Moon" reads retroactively like a child who has seen domestic violence in some way (or heard about it) but found a different outlet other than hurting physical beings. Such deflection often occurs with children of physically abusive parents; they learned somehow that physical violence is an "acceptable" reaction to certain negative feelings (acceptable in their household at least) but don't want to actually hurt a person, so they lash out at inanimate objects.
The more we see of the Diamonds, the more their dynamic becomes eerily familiar to those of us who come from homes with domestic abuse of various kinds, whether in the immediate family structure or in the generation previous. This feels more and more like a story of inherited trauma and abuse where the family in question was the ruling power of an entire society and thus the abuse shaped an entire culture. In short, the wider Homeworld culture had to stem from White Diamond, and Homeworld culture is one of oppression and abuse.
In this context, it's frightening to think what White may want to do to Steven or those who accompanied him. We've said before how Steven Universe will almost certainly lead to the reformation of gem culture to be more open and inclusive, and now that Steven's in a position of power on Homeworld, the beginning efforts will likely take place.