CRPropa 3.1 -- A low energy extension based on stochastic differential equations

TL;DR

This paper introduces CRPropa 3.1's new stochastic differential equation module for simulating low-energy cosmic ray propagation, enabling anisotropic diffusion modeling in complex magnetic fields with validated accuracy.

Contribution

The paper presents a novel low-energy propagation module in CRPropa 3.1 that uses stochastic differential equations to model anisotropic diffusion in cosmic ray transport.

Findings

Validated the accuracy of the new diffusion algorithm.

Demonstrated the impact of diffusion coefficient ratios on cosmic ray density.

Explored the influence of particle rigidity on diffusion behavior.

Abstract

The propagation of charged cosmic rays through the Galactic environment influences all aspects of the observation at Earth. Energy spectrum, composition and arrival directions are changed due to deflections in magnetic fields and interactions with the interstellar medium. Today the transport is simulated with different simulation methods either based on the solution of a transport equation (multi-particle picture) or a solution of an equation of motion (single-particle picture). We developed a new module for the publicly available propagation software CRPropa 3.1, where we implemented an algorithm to solve the transport equation using stochastic differential equations. This technique allows us to use a diffusion tensor which is anisotropic with respect to an arbitrary magnetic background field. The source code of CRPropa is written in C++ with python steering via SWIG which makes it easy to use and computationally fast. In this paper, we present the new low-energy propagation code together with validation procedures that are developed to proof the accuracy of the new implementation. Furthermore, we show first examples of the cosmic ray density evolution, which depends strongly on the ratio of the parallel $\kappa_\parallel$ and perpendicular $\kappa_\perp$ diffusion coefficients. This dependency is systematically examined as well the influence of the particle rigidity on the diffusion process.

Tables

Table 1. Table 1 : P-values for diffusion in homogeneous background field

	${10}^{12}$	$5.2\cdot {10}^{12}$	${10}^{13}$
X	0.64	0.95	0.06
Y	0.94	0.85	0.02
Z	0.23	0.58	0.30

Table 2. Table 2 : χ 2 superscript 𝜒 2 \chi^{2} /dof for diffusion in homogeneous background field

	${10}^{12}$	$5.2\cdot {10}^{12}$	${10}^{13}$
X	0.78	0.37	1.80
Y	0.38	0.48	2.22
Z	1.29	0.84	1.19

Table 3. Table 3 : Mean relative error for different simulation setups with fixed integration time h = 0.1 ​ pc / c ℎ 0.1 pc c h=0.1\,{\rm pc/c} .

	100	500	1000	2000
Z-boundary	0.023	0.045	0.048	0.049
Inner	-0.019	0.002	0.006	0.006
All	-0.011	0.011	0.014	0.015

Table 4. Table 4 : Mean relative error for different simulation setups with fixed number of snapshots N s ​ n ​ a ​ p = 500 subscript 𝑁 𝑠 𝑛 𝑎 𝑝 500 N_{snap}=500 .

	100	0.1-100	10	0.01-10	1	0.1
Z-boundary	0.554	0.188	0.172	0.12	0,045	0.008
Inner	0.113	0.033	0.036	0.019	0.002	-0.005
All	0.201	0.064	0.063	0.039	0.01	-0.003

Table 5. Table 5 : Mean deviation from field line in the JF12 field for different rigidities

$\rho \left[\mathrm{TV}\right]$	$10$	$100$	$1000$	${10}^4$	${10}^5$
$⟨\Delta R⟩\left[\mathrm{pc}\right]$	4.6	4.6	4.8	5.4	6.2

Table 6. Table 6 : Effective escape time scale for different diffusion parameter ϵ italic-ϵ \epsilon

$ϵ$	$1.0$	$0.1$	$0.01$	$0.001$
$T_{loss}\left[\mathrm{Mpc}/c\right]$	$19.2$	$23.3$	$29.4$	$32.6$

Table 7. Table 7 : Effective escape time scale for different rigidities ρ 𝜌 \rho and ϵ = 0.01 italic-ϵ 0.01 \epsilon=0.01

$\rho \left[\mathrm{TV}\right]$	$10$	$100$	$1000$	${10}^4$	${10}^5$
$T_{loss}\left[\mathrm{Mpc}/c\right]$	$29.4$	$14.1$	$6.9$	$3.6$	$2.0$

Equations

$$\frac{\partial n}{\partial t}+\vec{u}\cdot \nabla n=\nabla \cdot (\hat{\kappa }\nabla n)+\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{\kappa }_{pp}\frac{\partial n}{\partial p}\right)+\frac{1}{3}\left(\nabla \cdot \vec{u}\right)\frac{\partial n}{\partial \ln p}+S(\vec{x},p,t).$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}=a(x,t)+b(x,t)\xi (t),$$

$$⟨\xi (t+\tau )\xi (t)⟩$$

$$=\delta (\tau ).$$

$$x(t)=x(0)+{\int }_{t_0}^{t}a[x(s),s]\mathrm{d}s+{\int }_{t_0}^{t}b[x(s),s]\mathrm{d}W(s).$$

$$\mathrm{d}x(t)=a[x(t),t]\mathrm{d}t+b[x(t),t]\mathrm{d}W(t).$$

$$x(t_{i+1})=x(t_i)+a(x(t_i),t_i)\Delta t_i+b(x(t_i),t_i)\Delta W_i,$$

$$\frac{\partial n(\mathbf{x},t;\mathbf{y},t^{\prime })}{\partial t}=-\sum _i\frac{\partial }{\partial x_i}[A_i(\mathbf{x},t)n(\mathbf{x},t;\mathbf{y},t^{\prime })]+\frac{1}{2}\sum _{i,j}\frac{{\partial }^2}{\partial x_i\partial x_j}[B_{ij}(\mathbf{x},t)n(\mathbf{x},t;\mathbf{y},t^{\prime })],$$

$$\mathrm{d}{r}_{\nu }=A_{\nu }\mathrm{d}t+D_{\nu \mu }\mathrm{d}{\omega }^{\mu },$$

$$\frac{\partial n(\mathbf{x},t;\mathbf{y},t^{\prime })}{\partial t}=-\sum _i\frac{\partial }{\partial x_i}[A_i(\mathbf{x},t)n(\mathbf{x},t;\mathbf{y},t^{\prime })],\frac{\mathrm{d}x}{\mathrm{d}t}=A[\mathbf{x}(t)].$$

$$\mathrm{d}\vec{r}$$

$$\mathrm{d}q$$

$$D_{ij}={\delta }_{ij}\sqrt{2{\kappa }_{ij}},D_{qq}=\sqrt{2{\kappa }_{qq}}.$$

$${\vec{x}}_{n+1}$$

$$={\vec{x}}_n+\left(\sqrt{2{\kappa }_{\parallel }}{\eta }_{\parallel }{\vec{e}}_t+\sqrt{2{\kappa }_{\perp ,1}}{\eta }_{\perp ,1}{\vec{e}}_n+\sqrt{2{\kappa }_{\perp ,2}}{\eta }_{\perp ,2}{\vec{e}}_b\right)\sqrt{h},$$

$${\vec{e}}_t$$

$${\vec{e}}_n$$

$${\vec{e}}_b$$

$${\vec{r}}_{end}={\vec{r}}_{start}+{\int }_0^L\vec{B}/B\mathrm{d}s,$$

$$m=||{\vec{r}}_{end,4}-{\vec{r}}_{end,5}||\le \xi \cdot \mathrm{kpc},$$

$${\vec{r}}_{end}={\vec{r}}_{start}+\sum _{j=0}^{2^n-1}{\int }_{L_j}^{L_{j+1}}\vec{v}(s)\mathrm{d}s,$$

$$h_{next}$$

$$\frac{\partial n}{\partial t}$$

$$n(r_i,t)=\frac{1}{(4\pi {\kappa }_{ii}t{)}^{d/2}}\cdot \exp \left(-\frac{r_i^2}{4{\kappa }_{ii}t}\right).$$

$$\tilde{P}(R_i)$$

$${\chi }^2=\sum _{i=1}^n\frac{(O_i-E_i{)}^2}{E_i},$$

$$p({\chi }^2)={\int }_{{\chi }^2}^{\infty }{P}_k({\chi }^2)\mathrm{d}{\chi }^2,$$

$$-\nabla \cdot (\hat{\kappa }\nabla n(\vec{r}))$$

$$s(\vec{r})=\frac{4}{{\pi }^2}\left(\frac{2{\kappa }_{xx}}{4R^2}+\frac{{\kappa }_{zz}}{2H^2}\right)\cdot \cos \left(\frac{x\pi }{2R}\right)\cos \left(\frac{y\pi }{2R}\right)\cos \left(\frac{z\pi }{2H}\right),$$

$$n_{ana}(\vec{r})=\cos \left(\frac{x\pi }{2R}\right)\cos \left(\frac{y\pi }{2R}\right)\cos \left(\frac{z\pi }{2H}\right).$$

$$\frac{\partial n}{\partial t}=\nabla \cdot (\hat{\kappa }\nabla n)+s(\vec{r})\delta (t-t_i),$$

$$n_{sim}(\vec{r})=\sum _{i=1}^{N_{snap}}n(t_i)h_{i}w.$$

Full text

We developed a new module for the publicly available propagation software CRPropa 3.1, where we implemented an algorithm to solve the transport equation using stochastic differential equations. This technique allows us to use a diffusion tensor which is anisotropic with respect to an arbitrary magnetic background field. The source code of CRPropa is written in C++ with python steering via SWIG which makes it easy to use and computationally fast.

In this paper, we present the new low-energy propagation code together with validation procedures that are developed to proof the accuracy of the new implementation. Furthermore, we show first examples of the cosmic ray density evolution, which depends strongly on the ratio of the parallel $\kappa_{\parallel}$ and perpendicular $\kappa_{\perp}$ diffusion coefficients. This dependency is systematically examined as well the influence of the particle rigidity on the diffusion process.

1 Introduction

While the field of astroparticle physics has advanced significantly during the past decade coming closer to an identification of the sources, the origin of both Galactic and extragalactic cosmic rays is still not fully resolved. Modern experiments like IceCube [1, 2], Imaging Air Cherenkov Telescopes (H.E.S.S. [3], MAGIC [4], VERITAS [5]), the Pierre Auger Observatory (PAO) [6], the Telescope Array (TA) [7, 8] or AMS-02 [9] challenge theories and simulations with a high accuracy of the observed data. Especially the transition between the Galactic and extragalactic part of the energy spectrum is not fully understood. There are different simulation frameworks available for low-energy, Galactic and high-energy extragalactic propagation that meet a high standard already today [10, 11, 12]. One unified simulation framework, however, that can broadly and consistently describe all three major characteristics of the observations — namely chemical composition, energy spectrum and anisotropy level of the arrival directions — is not yet developed.

In particular, there are two different fundamental approaches for the modeling of cosmic ray (CR) transport: (1) The solution of the transport equation enables a approximate description in different environments, and thus has the major disadvantage that it is necessary to model the diffusion of the particles via a diffusion tensor, often reduced to a one dimensional diffusion coefficient. Here, also a variety of (semi-) analytic theories exist. These are e.g. the leaky box or more complex models (see e.g. [13, 14, 15, 16]). These models are able to describe the global shape of the spectrum but cannot fully reproduce the latest precision measurements. (2) The solution of the equation of motion, which is a very precise, but also a very CPU-time-consuming ansatz, only being usable at the highest energies, where it is feasible with respect to the computation times.

Transport Equation

On Galactic scales, the solution of a set of transport equations is commonly used, describing the flux for the different particle species and the change of their energy spectrum through processes like diffusion in space and momentum, advection and adiabatic cooling. A typical representation of the transport equation is given here:

[TABLE]

This is the so called Parker transport equation (augmented with a term describing momentum diffusion) which is a simplified version of Fokker-Planck transport equation. It is believed to be a good description of the particle transport problem in diffusive regimes (e.g. [10]). Here, $n$ is the particle density, $\vec{u}$ is the advection speed, $\hat{\kappa}$ is the spatial diffusion tensor, $p$ is absolute momentum, $\kappa_{pp}$ is the momentum diffusion coefficient used to describe re-acceleration and $S(\vec{x},p,t)$ describes the sources of cosmic rays. For the case of an entire nuclear network, as it is present in the cosmos, Equ. (1.1) can be extended to include catastrophic losses and gains from spallation (interaction with background matter) and nuclear decay. Furthermore, formulations of continuous losses from interactions with the magnetic field (e.g. synchrotron radiation) can be developed also for single species simulations.

Nowadays, there exist three major simulation tools for the Galactic cosmic ray transport: GALPROP [10], DRAGON [11] and PICARD [12]. GALPROP is probably the most widely used program but does not provide a method to use anisotropic, space dependent diffusion tensors, yet. Nonetheless, it is an excellent instrument to describe spallation and a variety of loss processes. DRAGON, which is based on the same principles as GALPROP, addresses the anisotropic diffusion problem among other things. It is known that the diffusive behavior along and perpendicular to the mean magnetic field line differ from each other (e.g. [17]). A more sophisticated numerical approach to solve the transport equation is PICARD. The biggest difference is that it uses multi-grids and provides a solver explicitly designed to yield the stationary solution of the transport equation.

Another technically different approach to solve the transport equation, commonly used for the propagation of cosmic rays inside the heliosphere [18, 19], makes use of the equivalence of parabolic partial differential equations and a corresponding set of stochastic differential equations, known from stochastic integrals and Ito calculus (see Section 2 and [20]). First attempts to apply this method also to the Galactic propagation of cosmic rays are made in [21, 22, 23, 24]. This concept is part of the methods introduced for Galactic propagation in this paper and will be discussed in detail in Section 2.

Equation of Motion

In CRPropa [25, 26] a single-particle approach is used for description of the transport of extra-galactic CRs. CRPropa 3.0 propagates single particles forward in time through the numerical integration of the equation of motion [26]. This opens up many new possibilities for the simulation. For example, no assumptions on the diffusion tensor have to be made, but arbitrary magnetic field configurations can be used. This makes the simulation more fundamental. In addition things like tracing of distinct particle paths and a stochastic treatment of interactions is possible. However, this method is not applicable to Galactic transport problems, because it is computationally too time consuming for simulations using particles with energies below a few $10\,{\rm PeV}$.

The energy range between the knee at  $10^{15}\,{\rm eV}$ [27] and the ankle at  $10^{18.5}\,{\rm eV}$ [6] in the cosmic ray spectrum is called transition region because there the source distribution shifts probably from Galactic to extragalactic sources (see e.g. [28] for a recent review). It is unclear where the influence of the extragalactic part of the spectrum begins to dominate. Furthermore, it would be interesting to decompose these two compositions in the knee-to-ankle region to learn more about possible sources with energies up to EeV energies. Up to now no open-access program exists which models the transition region as a combination of Galactic and extragalactic influences. Individual models have been presented describing a possible contribution of Galactic and extragalactic cosmic rays in this transition region [29, 30, 31, 32, 33]. In this paper, we introduce an extension of the public CRPropa code that enables the user to perform cosmic ray simulations at hundreds of PeV to ZeV energies with the single-particle propagation method, but also the modeling of TeV-PeV energies with the low energy add-on that we present in this paper using stochastic differential equations for the multi-particle propagation. This code is open-access111https://github.com/CRPropa/CRPropa3/releases and can be used by the entire community [26].

1.1 The Structure of CRPropa

In this section we explain the program structure of CRPropa 3.1 and how the new DiffusionSDE module is related to the other modules. CRPropa 3.1 is build of several independent modules which can be used to create a user defined simulation setup. This structure makes it easy to extend CRPropa with new modules, e.g. with new interaction modules.

The central data structure in CRPropa is the so called candidate. This object holds all necessary information about the (pseudo-) particle such as energy, direction, particle type, position and many more attributes. These attributes are altered and processed by the ModuleList during the propagation procedure. The ModuleList is basically a list containing all simulation modules chosen for a specific simulation. The user can load different kinds of modules into the ModuleList:

Boundary modules: Deactivate (stop simulation) the candidate if the boundary condition is met. 2. 2.

Source modules: Create new candidates and inject them into the simulation. 3. 3.

Interaction modules: Take interactions into account, e.g. proton-photon-interaction, radioactive decay, etc. 4. 4.

Deflection modules: Change the position of the candidate according to the equation of motion (PropgationCK) or transport equation (DiffusionSDE). 5. 5.

Observer modules: Detect particle if a certain condition is fulfilled. 6. 6.

Output modules: Save the simulation results to a text- or hdf-file.

Figure 1 shows the simulation process. The modules, except for the source module, are called successively and act on a candidate with each integration time step $h$ as long as the candidate is active.

The new module DiffusionSDE can replace the module PropagationCK as needed. In particular, we recommend to use the DiffusionSDE module for propagation with an energy below 100 PeV, while PropagationCK should be used for samples that are dominated by energies range above 100 PeV. While these numbers are a good estimate for the Galaxy they may differ significantly in other surroundings.

It cannot be generally answered which module is better. The two modules simulate two fundamentally different descriptions of the propagation process. Which module suits best depends on the specific simulation setup. In the case of weak magnetic fields with gyro radii of the order of the coherence length of the magnetic field a description by the transport equation is not possible. Also in the case of vanishing turbulent magnetic fields the particle transport is not diffusive and therefore can not be described by Equ. (1.1). In general, the distribution function $n$ has to be isotropic in momentum space to be described by the diffusion approximation.

On the other hand, the propagation of particles in very turbulent magnetic fields, where the gyro radius is smaller than the coherence length of the field, is computing time intensive. This may prevent the generation of sufficient statistics to draw significant conclusions.

In addition, we want to give a short overview of the different output formats and how they can be used to create observables which can be compared to measurements. CRPropa outputs can be generated in a tracking mode, where every (pseudo-) particle position during the simulation process is recorded. This is very useful if the concrete particle trajectory is of interest, but obviously yields huge amounts of data if it used for a full scale simulation with billions of candidates and hundred thousands of integration time steps for each candidate. In contrast to that, most other output modes save only those candidates which are detected by an ObserverModule. Here, the cosmic ray properties (position, momentum, etc.) at the origin and at detection are stored. With the DiffusionSDE module we developed also a new ObserverTimeEvolution module which takes snapshots of the total cosmic ray distribution for user defined time points.

Furthermore, not only cosmic rays but also neutral particles like $\gamma$ rays or neutrinos, which are created as secondaries during interactions, can be propagated and stored.

2 Transport Model and Stochastic Differential Equations

In order to extend CRPropa toward lower energies and by that to a dedicated Galactic propagation tool, we decided to use stochastic differential equations. This ansatz is most compatible with the single particle tracking used for the extra-galactic propagation.

Furthermore, it has some advantages over classical grid based methods. Since the calculation is not restricted to a grid, adaptive algorithms can be applied to vary the integration time step which depends on the local magnetic field geometry. This leads to a simulation as accurate as needed and as effective as possible. To take all possible interaction channels into account the grid based methods have to calculate the whole (DRAGON and GALPROP) or parts (PICARD) of the reaction network twice. That is not necessary following the stochastic approach, where all secondary particles are propagated as new particles.

The diffusion equation is solved in the local frame of the magnetic field line (see Section 2.3) which implies nearly no restriction on the allowed magnetic field set up. This means that no analytical description of the magnetic field is needed. In contrast to the approach in [22] the magnetic field direction is calculated on the fly and there is no need for derivatives of the vector field.

In times of huge CPU clusters the potential of parallelization is a crucial factor for a simulation program. The solution of the stochastic differential equation (SDE) is trivial to parallelize because the different phase space elements (hereafter called pseudo-particles) are independent of each other. The computation time for pure propagation of pseudo-particles scales linearly with the number of cores.

Furthermore, the simulation results can be reweigthed after the actual simulation is finished. The pseudo-particles can for example be weighted according to their energy at the source to change effectively the energy spectrum of the source distribution without the need of a new full size simulation. This makes it easy and computationally cheap to do parameter studies. It has to be mentioned that a change of the propagation parameters, such as diffusion coefficient or magnetic field configuration cannot be changed after the simulation.

2.1 Stochastic Differential Equations

In this section we briefly summarize the mathematical background of stochastic differential equations. Here, we only discuss the one-dimensional case which can in principle be generalized to more dimensions. Furthermore, it should be noted that the rigorous mathematical description is always just an approximation of a real physics problem, especially with respect to correlation times etc.

Stochastic differential equations have been examined since more than one hundred years ago. Starting with studies on the Brownian motion by Einstein and Langevin who derived the so called Langevin-equation independently:

[TABLE]

where $\xi(t)$ is a “rapidly [in time] fluctuating random term”, meaning $\xi(t)$ and $\xi(t^{\prime})$ are statistically independent for $t\neq t^{\prime}$ [20]. That means the random term is fully uncorrelated which can be described as:

[TABLE]

Furthermore, we require the time average to vanish ($\langle\xi(t)\rangle=0$) because any drift term can be described by $a(x,t)$. Equation (2.1) can consistently be interpreted by the corresponding integral equation (see e.g. [20]):

[TABLE]

Here, we used the fact that the integral of $\xi(t)$ is a Wiener process, $\int_{t_{0}}^{t}\xi(t^{\prime})\,{\rm d}t^{\prime}=W(t)$, which allows for the replacement of $\xi(s)\,{\rm d}s={\rm d}W(s)$. A Wiener Process is a Markov process and it can be interpreted as the solution of the diffusion equation with zero drift and a constant diffusion coefficient that is equal to 1. The sample paths are continuous but not differentiable.Equ. (2.3) is the prototype of a stochastic integral equation. The quantity $x(t)$ is determined by the initial condition $x(0)$ a deterministic term (first integral) and a stochastic term (second integral). The details of the existence and uniqueness of solution of these types of equations is beyond the scope of this paper and we refer the reader to the detailed studies in the textbook by Gardiner [20]. We just want to emphasize that the solution is not fully determined by the starting condition — as it is the case for ordinary differential equations — but also depends on the particular realization of the Wiener process. In addition it should be mentioned that the last term of Equ. (2.3) — the stochastic integral — can be interpreted in two ways, which differ in the choice where between the grid points $t_{i}$ and $t_{i-1}$ with a time ordering $t_{0}\leq t_{1}\leq\dots\leq t_{n}=t$, the stochastic part of the integral ($b[x(s),s)$] is evaluated. It can be shown that the stochastic differential equation is independent of the two, namely Itō and Stratonovich, interpretations of the stochastic integral [20].

If a stochastic quantity $x(t)$ is described by Equ. (2.3) for all $t$ and $t_{0}$ it obeys the Itō SDE given by:

[TABLE]

it can be concluded that the solution of the Itō SDE Equ. (2.4) is described by taking the limit $(t_{i+1}-t_{i})=\Delta t_{i}\rightarrow 0$ for the discretized version of Equ. (2.4):

[TABLE]

where $\Delta W_{i}=W(t_{i+1})-W(t_{i})$ and the grid points are time ordered $t_{0}\leq t_{1}\leq\dots\leq t_{n}=t$.

2.2 Connection between Fokker-Planck and Stochastic Differential Equations

The Fokker-Planck equation (FPE) is a parabolic partial differential equation and can be seen as a special case of the differential form of the Chapman-Kolmogorov equation with vanishing jump probability (e.g. [20]). In general, a FPE can be written as:

[TABLE]

where $A_{i}({\bf x},t)$ is the drift vector, $B_{ij}({\bf x},t)$ is the diffusion tensor and $n({\bf x},t;{\bf y},t^{\prime})$ is the density at place ${\bf x}$ and time $t$ depending on the density at place ${\bf y}$ and time $t^{\prime}$. Here, ${\bf x}$ and ${\bf y}$ are in principle higher dimensional phase space vectors. The propagation of the phase space elements is described by the corresponding stochastic differential equation (here for the four dimensional case with three spatial and one momentum dimension):

[TABLE]

where ${\rm d}t$ is the time increment, $r_{\nu}$ is a 4-dimensional vector $(\vec{r},||\vec{p}||)$ and ${\rm d}\omega^{\mu}=\sqrt{{\rm d}t}\;\eta^{\mu}$ stands for a 4-dimensional Wiener process with Gaussian noise. The square-root proportionality of ${\rm d}\omega$ to ${\rm d}t$ is easily connected to the general behavior of normal diffusion $\langle x^{2}\rangle/t\propto D$. $A_{\nu}$ is a transport vector, $D_{\nu\mu}$ is a 4x4 tensor representing the diffusion and $\eta^{\mu}$ is a 4-vector of independent normal distributed random variables with zero mean and unity variance.

The derivation of Equ. (2.7) from Equ. (2.6) is beyond the scope of this paper and the interested reader is referred to textbooks like e.g. [20]. The connection between the FPE in Equ. (2.6) and the SDE in Equ. (2.7) can be compared to Lioville’s equation and the ordinary differential equation which describes a deterministic motion:

[TABLE]

Indeed Liouville’s equation is only a special case of the Fokker-Planck equation for vanishing stochastic motion.

Since the spatial and momentum operators in Equ. (1.1) decouple for most Galactic applications, Equ. (2.7) can be split into two parts:

[TABLE]

where the absolute value of the momentum is denoted as $q=||\vec{p}||$, see e.g. [34] for a more detailed explanation.

The remaining problem is the calculation of the tensor $D_{r}$ and scalar $D_{q}$ in the SDE from the given spatial diffusion tensor $\hat{\kappa}_{r}$ and the given momentum diffusion scalar $/\kappa_{qq}$, respectively, according to Equ. (1.1). In general, one has to find the square root of the diffusion tensor $\hat{\kappa}$, which is the solution of $(\kappa+\kappa^{t})=DD^{\dagger}$. In our case this reduces to a simple uncoupled set of equations, because we neglect drift terms and solve the equation in the local frame of the magnetic field line. Hence, the diffusion tensor is always diagonal which leads two the following relations:

[TABLE]

Note that while the SDE method allows for the computation of ‘quasi-trajectories’ of particles and thereby, as mentioned above, makes it most compatible to the single particle tracking used for extragalactic cosmic rays, it nonetheless requires the knowledge of the diffusion tensor. This is opposite to so-called full-orbit simulations e.g. [35, 36] where the elements of the diffusion tensor are determined from the particle trajectories obtained from a direct solution of their equations of motion.

2.3 Numerical Implementation

In this section we will briefly explain the developed algorithm which calculates the next pseudo-particle position. First, we discuss the Euler-Mayurama scheme which can be used to integrate the SDE when the local orthonormal basis is known. In Section 2.3.2 the implemented adaptive method for the calculation of the local trihedron — orthonormal basis defined by the Frenet equation — is explained.

2.3.1 The Euler-Maruyama Scheme

To explain the diffusion code all terms except spatial diffusion are neglected in this section. Following Equ. (2.10) and Equ. (2.9) the so called Euler-Maruyama scheme (e.g. [37]) can be derived:

[TABLE]

where $h=t_{i+1}-t_{i}$ is the integration time step. Further, the orthonormal basis222The orthonormal basis is the local trihedron of the magnetic field line $\vec{B}(\vec{x}_{n})$ when it is interpreted as a three-dimensional curve. $\{\vec{e}_{t},\vec{e}_{n},\vec{e}_{b}\}$ is generally defined by the Frenet-Serret-equation (e.g. [38]) which can be rewritten for magnetic field lines (e.g. [19]) as:

[TABLE]

where $\vec{B}_{reg}$ is the regular or coherent background field vector and $k=|{\rm d^{2}}\vec{r}/{\rm d}s^{2}|$ is the curvature of the parametrized field line $\vec{r}(s)$.333In this paper the term “field line” always refers to the coherent background field vector $\vec{B}_{reg}$.. Generally, all three vectors of the local trihedron $\{\vec{e}_{t},\vec{e}_{n},\vec{e}_{b}\}$ have to be calculated on the fly at each time step. Although the two perpendicular diffusion coefficients differ in principle, we state for the Galactic scenario that they can be assumed equal $\kappa_{\perp,1}=\kappa_{\perp,2}$. Anisotropic perpendicular diffusion can, in general, either result from anisotropic (or non-axisymmetric) turbulence perpendicular to the mean magnetic field [39], which is not the case for the Galactic magnetic field (see the turbulent component of the JF12 field [40]). Otherwise the decoupling of the perpendicular directions is caused by a large curvature of the background field on scales of the mean scattering length of the particles. This is also not observed for the Galactic magnetic field. This is in contrast to simulations in the heliosphere, where the consideration of anisotrpic perpendicular diffusion can be crucial (see e.g. [41]).

2.3.2 Adaptive Calculation of the Local Trihedron

The tangential vector of the magnetic field line $\vec{e}_{t}$ is calculated via the so called Cash-Karp (CK) algorithm [42], an adaptive Runge-Kutta (RK) algorithm. With this algorithm we are able to adapt the time integration step to minimize the number of steps. The local truncation error for the field line integration can be set by the user. In this way the overall computation time is reduced without any losses in accuracy.

The tangential vector is approximated by the difference of two points on the magnetic field line $\vec{e}_{t}=\vec{r}_{end}-\vec{r}_{start}$. Here, the starting point is given by the old pseudo-particle position $\vec{r}_{start}=\vec{x}_{n}$ (see Equ. (2.12)) and the end position $\vec{r}_{end}$ is calculated with a field line integration:

[TABLE]

where $L=\sqrt{2\kappa_{\parallel}}\,\eta_{\parallel}\,\sqrt{h}$ is the arc length of the field line elements and has the length of the diffusion step in parallel direction. The CK algorithm uses a fourth and fifth order RK algorithm to calculate two solutions of the initial condition problem given in Equ. (2.14). The norm of the difference of the vectors is used as a measure for the local truncation error. The step is accepted if:

[TABLE]

where $\xi$ is the user set precision and $\vec{r}_{end,4}$ and $\vec{r}_{end,5}$ are the $4^{th}$- and $5^{th}$-order solution, respectively.

In the case of a rejected step integration length $L$ is bisected until the condition in Equ. (2.15) is fulfilled or the minimum step is reached. Following this procedure the end position vector $\vec{r}_{end}$, used for the calculation of the tangential vector $\vec{e}_{t}$, is derived by $2^{n}$ consecutive solutions of Equ. (2.14) if ($n-1$) step attempts are rejected:

[TABLE]

where $L_{j}=2^{-n}L\,j$ and $\vec{v}=\vec{B}/|\vec{B}|$ is the normalized magnetic field direction. This procedure is not as efficient as the conventional step adaption ($h_{next}=0.95\,h\cdot m^{-0.2}$ see [42]) but is necessary. If the step size is reduced as conventional, large values of $\eta_{\parallel}$ will be rejected more frequently as compared to small values. This would lead to an effective underestimation of the tails of the Wiener process ${\rm d}\omega$. Furthermore, this would not just be an inaccurate description of the diffusion process but is simply not an approximation of the diffusion. The error is hard to quantify since it depends on the global magnetic field structure, the diffusion tensor and also on the simulation set-up such as minimum/maximum step size and the precision. With the procedure of Equ. (2.16) we ensure that the particles are always propagated over the full drawn diffusion step but use an adequate number of intermediate steps.

The suggested next integration time $h_{next}$ is:

[TABLE]

When the current step is not decreased the next step should be increased. The next time ingration time step $H_{next}$ is increased by a factor 4, which leads on average to a doubled parallel step length $L$, due to the square root dependence of the step length from the integration time.

3 Validation of the Algorithm

The new code is validated in two ways. Firstly, we proof that our code reproduces the correct solution for a given diffusion coefficient. In doing so, we compare in Section 3.1 the simulation with simple analytical expectations. Further, we validate the adaptive field line integration in section Section 3.2.

All tests that are developed for this paper test a mathematical problem. That means they are not designed to describe the physical reality. So before the software is used for physical simulations it should be carefully examined if the necessary assumptions made for the tests do also hold (at least approximately) in the physical simulation.

3.1 Homogeneous Background Field

3.1.1 Green’s function

For the first test we simulate protons under the influence of diffusion in a homogeneous magnetic field $\vec{B}=B_{0}\cdot\vec{e}_{z}$. We use a diffusion tensor that is anisotropic with respect to the background field $\hat{\kappa}={\rm diag}(\kappa_{\perp},\kappa_{\perp},\kappa_{\parallel})$ with $\kappa_{\parallel}=10\kappa_{\perp}$. Here, the probability distribution function (PDF or density $n$) of a pseudo particle is known and given by the Green’s function of the following equation:

[TABLE]

which is solved by

[TABLE]

Here, $\kappa_{ii}$ is the diffusion coefficient, $d$ is the dimension ($d=1$ in this test case, because the three directions are decoupled) and $r_{i}$ is the $i^{\rm th}$ component of the spatial vector. In other words, we expect that the distributions of the three vector components ($x(t),y(t),z(t)$) of the pseudo-particles follow normal distributions with zero means and variances $\sigma_{i}=\sqrt{2\kappa_{ii}t}$ since that is the analytical solution as given in Equ. (3.2).

Figure 2 shows the density distributions for all three components at three different times. The difference between the parallel (z-) and the perpendicular (x-, y-) direction is clearly visible. Furthermore, it can be seen that the width of the distribution increases with increasing time and that the speed of this process depends on the diffusion coefficient. To emphasize this diffusive behavior we plot the absolute distance from the origin for all three components in combination with the analytical solution in Fig. 3. A difference between the analytical solution (line) and the simulated data (histogram) cannot be seen by eye. To quantify these results a statistical test is applied to the data.

Since we expect that the differences per bin between the simulated and expected distribution are normally distributed we use a $\chi^{2}$-test of the probability distribution function (PDF) of the particle range (norm of the end position vector components). So here, $P(r)$ is tested rather than $P(\vec{r})$. This PDF is described by:

[TABLE]

Here, $R_{i}=|r_{i}|$ is the absolute value of the vector components. The distance is sorted into bins while it is ensured that the minimal count per bin is greater than 5. This is done by a consecutive reduction of the number of bins. These bin counts are compared to the expected values derived from the integration of Equ. (3.3) over the bin edges. From that the $\chi^{2}$ value is calculated as:

[TABLE]

where $n$ denotes the number of bins and $O_{i}$ and $E_{i}$ the observed and expected count of events in the $i^{th}$ bin, respectively. Using the $\chi^{2}$ test result according to Equ. (3.4) we calculate the p-value:

[TABLE]

where $p$ is the p-value444Although the p-value is seen as deprecated by the American Statistical Association among others it is not critical to use it here. and $P_{k}(\chi^{2})$ is the $\chi^{2}$-distribution with $k$ degrees of freedom.

The Null hypothesis, that the sample is drawn from Equ. (3.3), is rejected if the p-value is not in the range ($0.005<p<0.995$). We choose this broad range of acceptance because we want to minimize the number of false rejections. Since problems with the algorithm, at least to our experience, lead to very large deviations we decided to use this robust acceptance range. This also allows us to use this test for the automated testing of CRPropa (see below). Here, a large number of statistical tests is applied and the test suite has to be rerun if one of them fails, so a very robust test setup is desired.

According to the p-values summarized in Table 1, we conclude that the algorithm reproduces the diffusive behavior correctly in the case of a homogeneous background field using an anisotropic diffusion tensor.

Unittest in CRPropa 3.1

A slightly different version of the tests described above is implemented in the unit test framework of CRPropa 3.1 which can be run within the installation process of the software. In this way it is ensured that the code is executed as expected. If all tests are passed, the user will be sure that no known software problems occur after the installation.

Here, two statistical tests are implemented. First, the Anderson-Darling test [43] is applied on the simulated data. This test can be compared to the commonly known Kolmogorov-Smirnov test but gives more weight to the tails of the distribution. The Anderson-Darling test is very sensitive to outliers of the data but cannot be used to test against a specific normal distribution. The Anderson-Darling turned out to be a very good test to find problems in the numerical implementation of the diffusion algorithm. Since is it sensitive to outliers the test detects even a small number of failed pseudo-particle propagations which was helpful during the development of the code.

If the Anderson-Darling test does not find any significant deviations from a normal distribution we can assume that the errors of the simulated data compared with the analytic results will follow also a normal distribution. This allows us to apply the $\chi^{2}$-test to verify that the width of the normal distribution does not deviate significantly from the expected one.

3.1.2 Stationary Test Solution

This test, as the one before, is not a physical model but rather a artificial problem specifically designed to test the algorithm. Therefore a comparison with real physical problems can be misleading.

In this section the test setup for a more complex simulation which is also used in the testing of the PICARD code [12] is explained. Here, the stationary solution of a simple diffusion equation:

[TABLE]

with $\hat{\kappa}={\rm diag}(\kappa_{xx},\kappa_{xx},\kappa_{zz})$ and a source term $s(\vec{r})$ is tested against the known analytical solution. With the boundary condition $n(x=\pm R,y=\pm R,z=\pm H)=0$ and the source term defined by

[TABLE]

where the solution is given by:

[TABLE]

For the test procedure we inject pseudo-particles using a simple rejection sampling method. In that way we make sure that the source distribution follows Equ. (3.7) ignoring the constant factor for now. Afterward they are propagated forward in time and the position of the particles is recorded at $N$ consecutive times $t_{1},t_{2},...,t_{N}$. The boundary condition is implemented simply by deactivating the candidates when they reach the boundary. This removal of the pseudo-particles from the simulation leads automatically to $n(x=\pm R,y=\pm R,z=\pm H)=0$ as required.

The time dependent particle distribution $n(t_{i})$ corresponds to the solution of the following equation:

[TABLE]

which differs from Equ. (3.6) since it is not a stationary equation. We assume that the particle distribution $n(t_{i})$ does not vary significantly in the time interval $h_{i}=t_{i+1}-t_{i}$ so we can calculate the stationary solution from our particle distribution as follows:

[TABLE]

Here, $w$ is a weight according to the source function, where $s(\vec{r})=w\cdot n_{ana}(\vec{r})$. In general, Equ. (3.10) is an infinite sum but the density $n(t)$ decreases with time $\lim_{t\rightarrow\infty}n(t)=0$. It follows that there exists a time $t=T_{\rm max}$ with $\int_{V}n(T_{\rm max})\,{\rm d}V\leq\beta N$ where only a mere fraction $\beta$ of the particle remain in the simulation volume. The maximum integration time $T_{\rm max}$ needs to be chosen such that $\beta$ is small enough for the aspired accuracy. In other words, since the time integration of Equ. (3.9) leads to Equ. (3.6), the discrete solution in Equ. (3.10) is an approximation for the solution of the original stationary transport equation.

This test is designed to answer several questions:

How many snapshots $N_{\rm snap}$ are needed to resolve the solution sufficiently? The number of snapshots $N_{snap}$ determines the number of summands in Equ. (3.10). 2. 2.

What is a reasonable size of the integration time step $h$? 3. 3.

What is the integration time $T_{\rm max}$ to reach the stationary solution?

In this test we implemented the parameters as $\kappa_{zz}=10^{24}\,\rm{m^{2}/s}$, $\kappa_{xx}=0.1\kappa_{zz}$, $R=H=0.5\,\rm{kpc}$, $T_{\rm max}=50\,\rm{kpc/c}$ and $N=10^{6}$. We used two different approaches: In the first one we fixed the integration time step $h=1\,{\rm pc/c}$ and varied the number of snapshots $N_{\rm snap}$. To answer the second question we fixed the number of snapshots $N_{\rm snap}=500$ and varied the integration time step $h$. In both setups we analyzed the integrated number density depending on the maximum integration time (see paragraph ‘Total particle number’). After that we examined the accuracy of the spatial density distribution $n_{sim}(\vec{r})$ for the latest time point $T_{max}=50\,\rm{kpc/c}$ (see paragraph ‘Spatial accuracy’).

Total Particle Number

Figure 4 and Fig. 5 each show the simulated integrated particle density $\int_{V}n_{sim}(\vec{r})\,{\rm d}V$ depending on the maximum integration time. In general both figures show that the approximated solution asymptotically increases to a value close to the expectation.

Figure 4 shows that the estimated density increases with an increasing number of snapshots. This is reasonable because a higher number of snapshots approximates the time evolution better. This gives effectively more weight to the early time points with higher densities as compared to the case with a low time resolution, where the first time snapshot is taken when already a significant number of particles is lost. The error for a sufficient number of snapshots $(N>500)$ is below one percent for $T_{\rm max}>50\,\rm{kpc/c}$. Furthermore, it is visible that the integrated particle number does not converge to a completely flat distribution, due to the fact that $\beta>0$. That means that not all pseudo-particles have time to leave the simulation volume before the end of the simulation time. A longer time integration can fix this error but is not computing-time efficient. Moreover, the density is slightly over-estimated. This is presumably caused by problems of the numerical implementation of the boundary conditions. To estimate the statistical error we repeated the simulation ten times and calculated the mean and standard deviation from these sets.

Figure 5 shows the results for different integration time steps $h$. Here, an increased accuracy (decreased time step) leads to a lower integrated particle number, where $h=0.1\,\rm{pc/c}$ gives an approximation smaller than the expected one. From that we conclude that an integration time step of $h=0.1\,\rm{pc/c}$ is small enough to fully resolve the boundary. In addition, we allowed the adaptive algorithm to choose the step size depending on the pseudo-particle position. In doing so, the algorithm becomes significantly faster than it would have been using the minimal allowed integration step only. On the other hand, it is a lot more accurate than the solution for the maximum allowed step.

Spatial Accuracy

After that the accuracy of the spatial density distribution $n_{sim}(\vec{r})$ is examined for $T_{max}=50\,\rm{kpc/c}$. In doing so, the data are binned in a three-dimensional histogram with ten bins in each direction to ensure sufficient statistics in each bin. Furthermore, the expected number of particles per bin is calculated using the bin edges and the analytic solution given in Equ. (3.8).

A first analysis showed that the error depends on the bin position. More precisely we can identify two groups of bins: The two outer layers in z-direction (orange bins in Fig. 6(a)) and the other bins (green bins in Fig. 6(b)). This anomaly of the errors is an effect of the anisotropic diffusion tensor. Since the mean spatial step is larger in z-direction than in the x-y-plane the boundary condition is resolved worse in these bins. In general, all simulations could reproduce the shape of the analytical solution. Figure 6(b) displays exemplary that the density distribution for $N_{snap}=500$ and $\Delta h=0.1\,\rm{pc/c}$ reproduces the expected shape of Equ. (3.8).

Table 3 shows the relative errors of setup 1 (fixed integration step) broken down for the different bin groups. The table gives also the overall mean per simulation setup. Evidently, the increase from 1000 to 2000 snapshots has nearly no effect on the results but it doubles the amount of data being processed. Consequently, a good guess for the required accuracy for the problem saves a lot of memory and processing time. The decreasing accuracy from 500 to 1000 snapshots can be explained by the fact that the approximated solution is not yet fully stationary, meaning the maximum integration time $T_{max}$ is not long enough.

Table 4 clearly exposes that a decreased integration time step $h$ decreases the relative error effectively. The good performance of the adaptive algorithm is confirmed also in the spatial analysis. The algorithm uses larger time steps in the middle of the simulation and decreases the integration time at the boundary. This results in a better accuracy than reached using the upper integration time limit at low computational costs. However, the accuracy of the lower time integration limit is not met.

A more detailed analysis, including full error resolution per bin and spread of the error in the different groups, is given in appendix A.

3.2 Field Line Integration

As it is shown in Section 3.1 the code reproduces the solution of the transport equation for a homogeneous magnetic field very well. The next step is to check whether the field line integration or derivation of the tangent vector $\vec{e}_{t}$ works. Therefore, we developed two different tests. The first problem has a simple analytic solution which enables us to calculate errors very carefully (see Section 3.2.1). The second one uses a realistic magnetic field as it may be chosen for a full simulation setup but does not provide easy analytic solutions (see Section 3.2.2).

As explained in Section 2.3.1 the algorithm is designed to adapt the integration time to follow the magnetic field line with the user defined precision. The basic assumption of this test is pure parallel diffusion, i. e. $\kappa_{\perp}=0$. Since magnetic field lines do not intersect, particles cannot leave their original field line. We make use of this fact in the following tests as we use the deviation from the field line as a measure for the accuracy of the algorithm.

3.2.1 Spiral Line

In this section we explain the test of the field line integration of the developed algorithm. A simple curve $\vec{r}_{\rm spiral}$ of the form:

[TABLE]

is used as the field line. The curve is parametrized in $z$ and the constant $s$ determines the windings per height. The form of the field line is shown in Fig. 7. Due to the varying curvature radius of the spiral Equ. (3.11) this curve is ideal to test the adaptive algorithm.

In doing so, we inject 10,000 particles at the origin $\vec{r}_{0}=\vec{0}$ into the simulation and propagate them forward for a time span of $T_{max}=100\,{\rm kpc/c}$. We used a minimum/maximum integration step of $h_{\rm min}=10^{-5}\,{\rm kpc/c}$ and $h_{\rm max}=1\,{\rm kpc/c}$, respectively. In different simulations the local truncation error $\xi$ is varied. The diffusion coefficient in parallel direction is fixed to $\kappa_{\parallel}\approx 6.4\cdot 10^{25}\,{\rm m^{2}/s}$ and the diffusion in perpendicular direction is neglected $\kappa_{\perp}=0$. The winding constant is fixed at $s=0.02\,{\rm kpc}$.

In addition to the position of the pseudo-particles we also tracked the arc length $L_{sim}$ of the particles trajectories. The simulated arc length $L_{sim}$ is calculated as the sum of the propagation steps considering the sign and not just the absolute value. From this variable we are able to reconstruct the true particle position. We calculate the arc length $L$ for an the end position $z_{max}=z(T_{max})$ and fixed start position $z_{0}=0$:

[TABLE]

Since Equ. (3.12) is not easily invertible we use a numerical solver to minimize $L_{sim}-L_{ana}(z)$. The best fitting $z_{ana}$ is used for the error estimation. Using $z_{ana}$ we can calculate the remaining two coordinates via Equ. (3.11). This deviation from the analytic position $\Delta_{1}=|\vec{r}_{sim}-\vec{r}_{ana}|$ is shown on the bottom left in Fig. 8.

Furthermore, we calculate the absolute distance from the field line as an alternative measure of the algorithm accuracy. In doing so, we simply minimize the distance between the simulated particle end position $\vec{r}_{sim}$ and the field line. This leads to the second error $\Delta_{2}=\min(|\vec{r}_{sim}-\vec{r}_{spiral}|)$ which is shown on the right of Fig. 8.

We repeated the test for precisions $\xi$ from $1$ to $10^{-9}$ and recorded the computation time for the total simulation555The simulation was done on a single core of the Intel®Core i7-6800 chip at $3.4\,{\rm GHz}$.. Figure 8 shows the results of this test. In the upper left the computation time is shown and it is clearly visible that there is no linear correlation between the precision and the computation time. This is due to the adaptive algorithm which uses small integration steps at low $|z|$ and increases the integration time in regions where the curvature radius is larger. In contrast, the deviation from the field line is nearly linearly correlated to the precision in the region where it is not bound by the maximum integration step (precision of $0.01-1$). The deviation from position $\Delta_{1}$ saturates not only at the low precisions but also at very high precisions. This is most likely due to the restricted resolution of the numerical minimizer.

This test clearly shows that the new adaptive algorithm is able to follow a field line with changing curvature radius with a user defined precision. Since the curvature radii of the Galactic magnetic field are in most cases much larger than in this test a precision of order $10^{-4}$ is most likely sufficient to receive sub-parsec errors from the propagation.

In general, it should be noted that such a large range of possible integration steps is not the most effective way to use an adaptive algorithm. A very high accuracy of order $10^{-7}$ or higher will probably lead to many (up to the order of 20 or more) consecutive step refinements. This can lead to millions of field line integrations for a single diffusion step (see Section 2.3.1 for details). The user should carefully choose a feasible range of steps and a reasonable precision. Good parameters are hard to predict but small benchmark simulations will likely help to find a useful set.

3.2.2 Galactic Magnetic Field

From the tests described in Section 3.2.1 it is known that the field line integration works. But what does that mean with respect to the Galactic environment? It is the aim of the proposed test here to evaluate whether the adaptive integration methods works on realistic Galactic field parameterizations.

The accuracy of the algorithm is supposed to be strongly correlated to the underlying magnetic field model. A simple continuously differentiable magnetic field is supposed to lead to a better result than a partly discontinuous field or one with poles. To get an estimate of the problems occurring from realistic Galactic magnetic fields, the magnetic field by Janson and Farrar (JF12-field) [40] was chosen for this test.

The magnetic field line is to be approximated with trajectory points along the magnetic field, again assuming $\kappa_{\perp}=0$. In doing so, magnetic field lines of $20\,{\rm kpc}$ length in $17\,{\rm pc}$ steps have been created for 100 randomly chosen emission points in the Galaxy.

In the following, a total number of $10^{5}$ particles are injected into the simulation. The start positions of the particles are randomly chosen among the 100 different sources locations for which ’best fit‘ field lines exist. The end position of the particles are recorded. Now the deviation from the field line is approximated by the minimum of the differences of field line points and end positions. In this way, 100 sets of minimum distances are created with a sample size of around 1000 each. This procedure is done for five different rigidities $\rho=(10^{13}-10^{17})\,{\rm V}$.

To analyze the data the mean distance $<\Delta R>$ is calculated for each rigidity. Table 5 shows that the error increases slightly with increasing rigidity. Higher rigidities imply increased mean spatial steps for the same time interval. The adaptive algorithm normally accounts for this problem by a reduction of the integration time. When the minimum integration time is reached the error cannot be reduced to the optimum which is visible here for the highest rigidities. The error does not converge to zero due to the finite resolution of the ideal magnetic field lines which leads to expected minimal error of $\Delta R_{\rm ideal}=4.6-4.7\,{\rm pc}$.

Although the current implementation of the JF12 field is not completely continuously differentiable the adaptive algorithm works well.

4 Cosmic Ray Density Evolution

There are many possible scenarios to use the new code, starting with the simulation of nearby known supernova remnants (SNRs) or a backtracking of particles, originating in the solar system to identify possible source regions inside the Galaxy. We decided to perform a different approach and simulate the diffusion of cosmic rays in a global sense. In order to get an estimate of the overall cosmic ray distribution we use a continuous source distribution following the SNR distribution given in [44]. Furthermore, the realistic JF12 field, explained in Section 3.2.2, is used as the regular background field for the diffusion. Different values for the diffusion parameter $\epsilon$ with $\kappa_{\perp}=\epsilon\kappa_{\parallel}$ are tested, as well as the influence of the rigidity on the density evolution.

Since we fix the energy for each simulation we do not have to include any process that might change the energy of the pseudo-particles. For a realistic, full-scale simulation of the CRs in the Galaxy this should of course be included. These processes are for example energy losses by interaction but also diffusive re-acceleration, which could in principle be implemented as a diffusion process in momentum also using stochastic differential equations.

Before we explain the results in Section 4.3 the source distribution and the simulation setup are explained in Section 4.1 and in Section 4.2, respectively.

4.1 Source distribution

We use the same supernova distribution as in [44] where the radial distribution is originally given by Case and Bhattacharya [45]. Although this distribution just gives a rough estimate it should reflect the general features of the SNR distribution inside our Galaxy well enough. The radial distribution is described by:

[TABLE]

where $\beta=3.53$ and $R_{0}=8.5\,{\rm kpc}$ are taken from [44] and $A_{R}=\beta^{4}\exp(-\beta)/(12\pi)$ is the normalization constant. The distribution in z-direction is defined by:

[TABLE]

where $z_{g}=0.3\,{\rm kpc}$ is chosen as the scale height of the distribution and $A_{Z}=1$ is once more a normalization constant.

To generate the pseudo-particle start position we use a rejection sampling for $r$ and $z$ from Equ. (4.1) and Equ. (4.2), respectively. The third coordinate in cylinder coordinates, the angle $\phi$, is drawn from a uniform distribution with $\phi\in(0,2\pi]$. After that the coordinates are transformed into cartesian coordinates $x=r\cos(\phi)$, $y=r\sin(\phi)$ and $z=z$. Our approach is different from the one in [44] since we do not define a given number of SNR positions drawn from the distribution but every particle has its own, unique origin. The sampling technique implemented for the source generation is very simple but fast enough. We can generate one million source positions in a few seconds, meaning that the generation of sources is not a significant part of the overall simulation time.

Figure 9 shows the distribution of SNRs used in the simulation. It is visible that the density has a minimum at $r_{min}=0$ and and maximum around $r_{max}=2R_{0}/\beta\approx 4.8\,{\rm kpc}$. After that the density decreases exponentially and is cut off at the boundary of the simulation volume at $r>20\,{\rm kpc}$ (see Section 4.2). The rotation symmetry is clearly visible which should be kept in mind to interpret the results in Section 4.3.

4.2 Simulation Setup

In this section the other simulation parameters, beside the source distribution, are discussed in detail.

As mentioned above, the regular magnetic field component of the JF12 field [40] is used as the background field. Due to the fact that we only use the regular component we have to slightly adapt the field. In the original implementation of the field the regular field component of the field vanishes for $|\vec{r}|\leq 1\,{\rm kpc}$, where the origin $\vec{r}_{0}=0$ is at the Galactic center (meaning $B_{reg}=0$ in this region). Since this zero-field sphere has no physical reason but is due to a lack of precise observational data we continue the field to $\vec{r}_{0}=\vec{0}$. In doing so, we neglect the hard cut off at $|\vec{r}|=1\,{\rm kpc}$ which is possible since the analytic form of the field is continuous up to $\vec{r}_{0}$.

This simulation is primarily designed to study the influence of the diffusion parameter $\epsilon$ and rigidity $\rho$ on the cosmic ray distribution. The diffusion tensor $\hat{\kappa}$ is diagonal in the local frame of the magnetic field line, $\hat{\kappa}={\rm diag}(\kappa_{\perp},\kappa_{\perp},\kappa_{\parallel})$ with the parallel diffusion coefficient defined as in [10]:

[TABLE]

where $\kappa_{0}=6.1\cdot 10^{24}\,{\rm m^{2}/s}$, $\alpha=0.3$ and $A(\epsilon)=1.02/(1+2\epsilon)$ is a normalization constant. The normalization $A$ is chosen to keep the trace of the diffusion tensor constant for different values of the diffusion parameter $\epsilon$. This is necessary since the value of the trace has a huge impact on the loss time scale. Using this normalization we exclude the influence of a varying trace on the time scale and are able to analyze the influence of the diffusion parameter $\epsilon$ or rigidity $\rho$.

In every simulation setup we calculate the trajectories for about 15 million pseudo particle forward in time until the maximum integration time $T_{max}=100\,{\rm Mpc/c}$ or the escaping from the simulation volume at $|\vec{r}_{max}|=20\,{\rm kpc}$. We use adaptive steps with $h_{min}=0.1\,{\rm pc/c}$, $h_{max}=1.0\,{\rm kpc/c}$ and precision $\xi=10^{-5}$. In that way we ensure the desired accuracy for the averaging process (see below). The particle position is recorded at 1000 points in time $t_{n}=n\Delta t$ with $\Delta t=100\,{\rm kpc/c}$ and $n=1...1000$. For the first analysis we use protons with a rigidity of $\rho=10\,{\rm TV}$ and for the second analysis we iteratively increase the proton rigidity by an order of magnitude up to $\rho=100\,{\rm PV}$.

The simulation is performed on a cluster with 256 cores and takes about 3200 CPU hours per simulation set. The raw output files take a total memory storage of up to 300 GB per simulation set, depending on the escape time scale. To handle this huge amount of data we bin the pseudo-particle positions into a three-dimensional histogram with $1\,{\rm kpc}$ bin edges for each time step $t_{n}$ and simulation setup. This concurrently ensures that the particle density is averaged over a sufficiently large number of pseudo-particles.

4.3 Results

The first results show that the diffusion along the magnetic field lines affect the cosmic ray density evolution significantly. We can show that even a nearly unstructured source distribution is transformed into a complicated density distribution following the magnetic field structure of the Galaxy. This is in contrast to earlier simulations where the source distribution itself is already shaped in a spiral structure e.g. [10, 12]. We discuss the results for varying diffusion parameter $\epsilon$ (Section 4.3.1) and rigidity $\rho$ (Section 4.3.2) separately.

4.3.1 Diffusion Parameter

Figure 10 shows different face-on views (anti-parallel to the z-axis) of the Galactic disc. In order to highlight the Galactic structure we display only the particles near the Galactic plane with height $|z|\leq 1\,{\rm kpc}$. The difference in the diffusion parameter $\epsilon$ is clearly present. As expected, a decreased perpendicular diffusion coefficient emphasizes the structure of the magnetic field. In the case of very strong perpendicular diffusion $\kappa_{\perp}=\kappa_{\parallel}$ almost any structure, besides a fall off in the direction of the boundary, can be found in the density distribution.

Another interesting fact is the difference in the cosmic ray density at the Galactic center. A strong perpendicular diffusion $\epsilon\geq 0.1$ leads to an increase of the Cosmic ray density. Finally, the maximum density region is shifted into the Galactic center. In contrast, strong parallel diffusion prevents the pseudo-particles from diffusion into the central region.

In addition to the analysis of the Galactic disc, we examined the cosmic ray density in the halo as well. In doing so, we combined in Fig. 11 all particles to produce edge-on representations of the Galaxy. As in the face-on projection the unstructured state of the cosmic ray density for (very) strong perpendicular diffusion is clear. The general trend of higher escape rates in the case of increasing perpendicular diffusion is apparent and quantified in Fig. 12.

Furthermore, especially for strong parallel diffusion, a remarkable difference between the northern and southern hemisphere can be noticed. In the North the cosmic rays are concentrated in the center for lower Galactic latitudes. Perpendicular outflows are produced as clearly visible in the figures by the halo component of the magnetic field. These outflows have shapes that resemble the x-shape structure of the magnetic field that is detected for external galaxies (e.g. [46, 47]) and their winds (e.g. [48]). In this paper, we do not perform a dedicated analysis in this direction and cannot draw any physics conclusions from these first findings. Detailed studies on the role of diffusion and advection for galactic outflows are planned in the future, but go beyond the scope of this paper. A similar structured outflow is not observed for the southern part of the Galaxy. This difference is an effect of the different magnetic field structure of the JF12 field in both hemispheres because the radial transition width of the toroidal field cannot be constrained by data for the southern hemisphere (see [40, 49] for details). This broken symmetry is not necessarily a real effect: the lack of observations for the halo field in the southern hemisphere leads to an approach to use the most simple magnetic field representation. Thus, even in the south, there might exist an x-shape-like structure, but it cannot be resolved by data yet at this point.

The escape time scale of the different diffusion configurations was mentioned above and is quantified in Fig. 12. This is not surprising, although the different setups have equal overall diffusion strength, but describe a completely different morphology. This difference is mainly caused by the spiral magnetic field component which is strongly aligned in the Galactic plane. Since large perpendicular diffusion coefficients $\kappa_{\perp}\approx\kappa_{\parallel}$ allow a fast escape into the halo these configuration have shorter loss times. In contrast, a dominant parallel diffusion component $\kappa_{\parallel}\gg\kappa_{\perp}$ binds the pseudo-particles strongly inside the Galactic disc preventing a fast escape.

We define the effective escape time scale $T_{loss}$ by

[TABLE]

as the time scale on which the total particle density decreases by $1/e$. The different time scale are given in Table 6.

4.3.2 Rigidity dependence

In this section the influence of the particle rigidity on the evolution of the cosmic ray density is discussed. We fix the diffusion parameter $\epsilon=0.01$ and vary the rigidity on a logarithmic scale from $\rho=10-10^{5}\,{\rm TV}$. The analysis is the same as described in Section 4.3.1. High rigidities lead to a very fast particle loss, so we decided to show the time evolution on a logarithmic scale rather than on a linear scale (as in Section 4.3.1).

Figure 13 shows the density distribution in the Galactic disc at $t=(0.1,1,10,100)\,{\rm Mpc/c}$ for different rigidities. The very left column may help to compare the two analyses as the data for $\rho=10\,{\rm TV}$ was already used in the previous analysis. It is visible that the simulation for different rigidities do not differ much from each other when the general topology of the density distribution is compared. But it is obvious that an increase in the particle rigidity leads to an acceleration of the escape process. The different simulations are very similar to each other but on very different time scales.

The effect of different time scales becomes even more prominent in the edge-on view shown in Fig. 14. At the highest rigidity the particles fill up nearly the whole halo in just a few hundred thousand years. The early time evolution reveals a feature of the cosmic ray transport that was not visible in the first analysis. The x-shape like structure which is observed for the northern hemisphere at later time is also visible in the southern hemisphere but on much shorter time scale. After a couple million years, depending on the rigidity, the diffusion in the toroidal magnetic field washes the x-shape structures out.

Figure 15 emphasizes the huge differences in the escape time scales. for rigidities above $\rho=1\,{\rm PV}$ it takes less than a hundred million years to lose nearly all particles to the extra-galactic medium.

In Table 7 the escape time scales are listed for a quantified analysis. Comparing the time scales we obtain that

[TABLE]

which reflects the rigidity dependence of the diffusion coefficient given in Equ. (4.3).

5 Summary and Outlook

In this section we emphasize once more the advantages of the new algorithm and recapitulate what is already possible using SDEs for Galactic propagation. In addition, we give a short outlook into the further development of simulation tasks and possible improvements of the status quo.

In principle the diffusion approach as it is implemented here is also able to simulate CRs with energies below $E\leq 10\,{\mathrm{T}eV}$. The problem here is that CRPropa uses approximations for very high energies when it comes to interactions processes. For example secondary particles are ejected parallel to the parent nucleus. This makes a low energy approach not impossible but surely the interaction modules have to revised.

5.1 Summary

In this paper we show that stochastic differential equations provide an efficient mathematical tool to describe the cosmic ray transport in the diffusive regime. This method is not a completely new approach but is already established in a few propagation codes used for the simulation in the heliosphere (e.g. [19, 34]). For Galactic propagation only first examples exist (e.g. [22, 23]). But in contrast to the approaches mentioned above we do not solve the transport equation in the global or laboratory frame. As explained in detail in Section 2.3 we take advantage of the fact that the diffusion tensor is diagonal in the frame of the coherent magnetic background field. This allows us to decouple the diffusion process for the three distinct directions of the orthormal basis defined by the comoving trihedron $\{\vec{e}_{t},\vec{e}_{n},\vec{e}_{b}\}$ (see Equ. (2.13)). In doing so, the diffusion tensor is completely defined by two parameters, the parallel diffusion coefficient $\kappa_{\parallel}$ and the diffusion ratio $\epsilon$.

We prove that the adaptive field line integration works in simple fields (Section 3.2.1) as well as realistic field descriptions (Section 3.2.1). Furthermore, a method to yield stationary solutions of the transport equation by averaging the time-dependent solutions is given and explained in Section 3.1.2.

In Section 4 first basic simulation examples show the potential of this new simulation software. Already this very simple source distribution leads to interesting insights of the time evolution of the cosmic ray density in the Galactic halo. These need to be quantified and analyzed in future work.

We conclude that this software may help to develop a code to describe the transition region consistent in a single propagation framework. In this way systematic differences between the treatment of Galactic and extra-galactic cosmic rays can be minimized.

Finally, we want to emphasize that this propagation tool is restricted to Galactic propagation but can be used in any environment where the diffusion approximation holds. In principle, this technique may also be used for the simulation of reacceleration by the implementation of momentum diffusion.

5.2 Outlook

In this section we will give a brief outlook on the future development and possible applications of the described tool.

In the case of our own Galaxy we will examine the cosmic ray density evolution further. In doing so, we will use different source scenarios, e.g. applying discrete SNRs and not a continuous source distribution or more sophisticated source models taking the spiral structure into account (e.g. the distribution of isolated radio pulsars [50]).

In the future this new program will contribute to the solution of the gradient problem of the observed gamma-rays (e.g. [51]) with the possibility of combining dedicated propagation tools with precise interaction models. The density maps shown in Section 4 are first step into the modeling of the cosmic ray gradients. The implemented anisotropic diffusion tensor will also help to understand the anisotropy of cosmic rays better.

From the technical point of view, several major extensions of the existing software are planned. The new propagation software has some advantages over the grid-based methods as mentioned above. But up to now the software is not able to take all kinds of interactions, like proton-proton or proton-nucleon inelastic scattering processes into account. Only proton-photon interaction are already implemented very efficiently in the framework. This means in effect, a complete description of the transport of cosmic rays in a dense environment like our Galaxy is not yet possible. The implementation of the missing interaction processes is one of the next steps to improve the code. In doing so, we will adapt the methods already used and tested in CRPropa for proton-photon scattering. Here, we will use Monte-Carlo methods to generate secondary particles depending on the specific cross sections (see e.g. [52, 53]) and inject them as new candidates into the simulation chain. Furthermore, density maps of the Galactic mass distribution have to implemented for a proper, location-dependent simulation of the spallation processes. Beside the implementation of the scattering processes three dimensional density maps of the target material (e.g. HI and HII maps from [54]) are crucial to derive realistic maps of neutral secondaries, such as gamma rays and neutrinos. We want to emphasise that implementation of interactions is relevant to model the fluxes of secondary particles like lithium, beryllium, boron, nitrogen, electrons, positrons and neutrinos. On the other hand the flux of primaries, such as protons, will not be affected much by interactions because the interaction probability is much smaller than unity for the majority of the CR population. An exception would be the combined case of strong parallel diffusion ($\epsilon\leq 0.01$) and low rigidities ($\rho\leq 100\,{\rm TV}$) where primary fluxes can be significantly influenced by interactions.

To take Galactic winds or in general advection processes into account an additional term of the transport equation has to be factored in. The advection is described by the linear term $A$ in Equ. (2.7) which means that the Euler-Maruyama Scheme in Equ. (2.12) is modified by a further term proportional to the integration time step $h$.

Up to now the diffusion tensor is locally independent in the frame of the magnetic field line. This means that the diffusion tensor, transformed in the lab frame, has constant eigenvalues. Recent studies by Snodin and others have shown that the calculation of the diffusion tensor depending on the local magnetic field structure and the particle’s gyro radius is possible [55]. They used test particle simulation to derive the diffusion tensor for a homogeneous background field and different models for the random magnetic field component. They give analytic formulas for the parallel and the perpendicular diffusion coefficient which might serve as a good starting point for a spatially varying implementation of the diffusion tensor in CRPropa 3.1.

The technically most challenging problem we try to solve is the development of a dynamic switching between the two different propagation models. As mentioned in Section 1.1 the definition of parameters for the change between the propagation modes is non trivial. We will use machine learning algorithms to find the best simulation procedure depending on the current candidate properties, like particle type, position and rigidity. First studies may for example compare the escape times in equal environments using the different modules. The long term goal is to develop a software which is able to dynamically choose the best propagation module depending on user set criteria. These criteria should take limited computing resources as well as the specific physical simulation into account.

Appendix A Spatial error resolution stationary solution

Acknowledgments

The authors would like to thank the whole CRPropa development team, especially G. M ller, A. Dundovic and D. Walz, for their persistent help with the implementation of the algorithm. We also appreciate discussions with A. Kopp and M.S. Potgieter. Some of this work was made possible by the following python packages: NumPy and SciPy [56], Pandas [57], matplotlib [58] and IPython [59]. We would further like to acknowledge support from the RAPP Center (Ruhr Astroparticle and Plasmaphysics Center sponsored by the MERCUR project St-2014-0040) and from the Research Department of Plasmas with Complex Interactions (Bochum). GS is supported by the DFG through collaborative research centre SFB 676, by the Helmholtz Alliance for Astroparticle Physics (HAP) funded by the Initiative and Networking Fund of the Helmholtz Association and by the Bundesministerium f r Bildung und Forschung (BMBF) through grant 05A14GU1.