{'fraction_non_alphanumeric': 0.03521842585533763, 'fraction_numerical': 0.01732621803312773, 'mean_word_length': 4.458701134930643, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmic-ray propagation with DRAGON2: I. numerical solver and astrophysical ingredients

We present version 2 of the DRAGON code designed for computing realistic predictions of the CR densities in the Galaxy. The code numerically solves the interstellar CR transport equation (including inhomogeneous and anisotropic diffusion, either in space and momentum, advective transport and energy losses), under realistic conditions. The new version includes an updated numerical solver and several models for the astrophysical ingredients involved in the transport equation. Improvements in the accuracy of the numerical solution are proved against analytical solutions and in reference diffusion scenarios. The novel features implemented in the code allow to simulate the diverse scenarios proposed to reproduce the most recent measurements of local and diffuse CR fluxes, going beyond the limitations of the homogeneous galactic transport paradigm. To this end, several applications using DRAGON2 are presented as well. The new version facilitates the users to include their own physical models by means of a modular C++ structure.

detected fluxes. DM indirect searches are strictly affected by uncertainties on the propagation of CRs in the Galaxy. Therefore, reducing and quantifying these astrophysical uncertainties is crucial in order to be able to disentangle a possible DM signature from astrophysical backgrounds and to determine the experimental sensitivity to DM properties.
In the wake of the plethora of novel experimental achievements described above, a parallel effort on the CR-transport modelling in the Galaxy is needed [25].
The DRAGON project has been pursued in order to meet the urgent demand to model CR propagation under the most realistic and general conditions.
Within this context, we introduce here the code version 2 (DRAGON2) as a general tool to simulate all relevant processes regarding CR transport from very low (∼ 10 MeV and below) to extremely high (∼ 1 PeV) energies. 1 In particular, it computes the solution of the diffusion-advection-loss equation describing CR transport for most CR species, from heavier ones down to protons, antiprotons, and leptons, both of astrophysical and exotic origin (i.e., coming from DM annihilations/decays). The transport equation features fully position-and energy-dependent transport coefficients (spatial and momentum diffusion, energy losses and advection) in both spatial two-dimensional (assuming cylindrical symmetry) and three-dimensional mode.
DRAGON2 allows a detailed study of both small-scale and large-scale structures (e.g., the spiral structure of the Galaxy) in steady-state and transient mode, refining the spatial resolution on the regions of interest (e.g., local bubble, GC, or Galactic Plane).
In this paper we introduce the code having in mind the specific case of CR propagation in the Galaxy; however, the code is written in a general way and can be easily used in many other different contexts and on different scales (e.g., CR transport in Galaxy clusters or in a star forming region).
Moreover, the new modular structure of the code make it possible for the user to implement additional spatial distributions for all relevant astrophysical quantities in a straightforward way.
Old versions of DRAGON [26] have been used in several contexts. For example, to provide a solution to the CR gradient and isotropy problem in terms of inhomogeneous diffusion [27]; to compute CR electron and positron spectra in the presence of a spiral arm structure for sources [28]; to model CR antiprotons as CR secondaries and from DM annihilations [29,30]; to study the synchrotron emission from galactic leptons [31]; and to reproduce γ-ray and neutrino diffuse emissions above the TeV [32]. Most of these original results were obtained thanks to the innovative features already present in the first version of the code: these features are included also in the new version of the code and, where possible, extended to more general cases (see the discussion in Section 5).
With respect to the previous versions, the new code DRAGON2 was also largely reworked in order to optimally profit from modern programming design and computing techniques. In this paper, we provide a detailed description of the transport equation solver (see Section 3) and provide in the appendixes details about different models for the relevant astrophysical ingredients (e.g., source and gas distribution, magnetic field models, spiral arm patterns) adopted (see Appendix C). In Section 4, we also provide a comprehensive set of numerical tests to assess the code performances, and to study the accuracy of the solution and the time needed to reach convergence in different conditions (e.g., different grid sizes, constant and variable time step). The main new features for galactic propagation are demonstrated in Section 5 in a few example applications.
The first large project of this kind, GALPROP 2 , is a widely used code in the community [33][34][35]. GALPROP is designed to make predictions of direct CR measurements as well as gamma rays and synchrotron radiation consistently. It includes realistic models for nuclear spallation processes [36][37][38][39][40] and energy losses, but basic assumptions for the CR transport 3 . Semi-analytical solutions of the propagation equation are implemented in the USINE code developed since 2010 [41]. Taking advantage of much faster computation methods than numerical models, the semianalytical approach allows for an efficient scan of a wider transport parameter space [42][43][44].
Recently, the PICARD numerical code have been developed [45,46]. PICARD is fully 3D in concept and implements modern numerical techniques for the numerical solver, handling high resolutions with reasonable computer resources.
DRAGON2 is part of a complete suite of numerical tools designed to cover most of the relevant processes involving Galactic CRs and their secondary products over a very wide energy range. With the help of these tools -in particular the HeSky 4 package -it is possible to compute spectra and sky-maps of radiation emitted by CRs interactions in a huge energy range, from the synchrotron radio waves up to the PeV neutrinos. On the low-energy side, the solar modulation can be treated either with auxiliary analytical routines implementing the force-field approximation [47,48], or with the HelioProp numerical code featuring a detailed model of CR charge-dependent interaction with the Heliosphere, including diffusion, advection and energy losses due to the solar wind [49]. This paper does not contain a description of spallation processes and of off-diagonal anisotropic diffusion, which will be covered in forthcoming publications and in the evolving DRAGON manual (see www.dragonproject.org).

Transport of CRs in the Galaxy
DRAGON2 features all relevant processes for CR transport from Galactic acceleration sites to Earth: in particular, spatial and momentum diffusion, energy losses, advection, nuclear spallations and decays.
The combination of all these processes can be described by the following equation [50,51]: where N i ( r, p) is the density per total momentum p of the CR species i, D pp ( r, p) is the momentum diffusion coefficient, Q( r, p) describes the distribution and the energy spectra of sources, v w ( r) is the Galactic wind velocity responsible for CR advection,ṗ( r, p) accounts for the momentum losses. The timescale for radioactive decay at rest is given by τ i , while σ i is the spallation cross-section with the interstellar gas. In this paper we do not consider these latter nuclear processes, and we postpone a detailed description to a forthcoming publication. The CR macroscopic current J( r, p) is determined by the spatial diffusion tensor D ij , as J i = −D ij ∇ j N . These quantities can be either inferred from independent observations (e.g. the gas distribution, the magnetic field entering the loss term) or fitted to the data (e.g. the diffusion coefficient, the Galactic wind velocity). For all of them, different parameterizations are provided in literature and can be used to estimate the systematic uncertainty affecting the corresponding process. We therefore implement in DRAGON2 several options for the relevant transport quantities, as extensively described in Appendix C; in most cases, the quantities are position-dependent.
As discussed in the Introduction, one of the main novelty of our code with respect to other existing codes is the possibility to implement inhomogeneous transport 5 (e.g., advection, momentum and spatial diffusion).
In particular, assuming diffusion as inhomogeneous and anisotropic has a very natural motivation. In fact, the presence of a large scale Galactic magnetic field (GMF) clearly breaks isotropy and introduces a preferred direction, so that charged-particle diffusion should be expressed in terms of a diffusion tensor with components given by: where b is a unit vector along the mean (large scale) GMF. With this choice of versors, D and D ⊥ are the components of the diffusion tensor parallel and perpendicular to the mean magnetic field and describe diffusion due to small-scale turbulent fluctuations. The coefficient D A gauges the anti-symmetric component of the diffusion tensor: It is usually identified as the drift coefficient since it describes a macroscopic drift orthogonal to both b and the gradient of the CR density, ∇N [52,53]. In this paper we always assume D A = 0 since the associated drifts are negligible up to ∼PeV energies as shown, e.g., in [54]. Although the physics behind CR diffusion is far from being understood (see e.g. [55] for a comprehensive review), some basic aspects may however be clarified starting from the weakturbulence approximation where magnetic perturbations are well-developed in k-space and small compared with the regular background component. Under this assumption it is possible to treat analytically the problem of resonant CR interactions with the random-phase MHD wavemodes. This framework is known as quasi-linear theory (QLT) [56,57]. The classical result for QLT gives that diffusion coefficients are described by a power-law in rigidity with different slopes for the parallel and perpendicular components (see also [58]). Moreover, these coefficients are spatially inhomogeneous since they are determined by local properties of the turbulent and regular fields. In this perspective, for the diffusion coefficients D and D ⊥ we adopt several phenomenological parameterizations as proposed in recent works based on local fluxes and gamma-ray data (see Appendix C.8).
DRAGON2 can work either in a (2 + 1)-dimensional (2D) or in a (3 + 1)-dimensional (3D) configuration. In the 2D case we use cylindrical coordinates defined by the radial distance r and the height form the Galactic disk z and we assume azimuthally symmetry. For the 3D case we consider Cartesian coordinates x, y, z. The quantities defined as function of cylindrical coordinates are consistently mapped in Cartesian coordinates by the relation r = x 2 + y 2 .
In the next Sections, we will specify the transport equation in these two configurations.

2D in cylindrical coordinates and azimuthal symmetry
In 2D mode, D plays no role and the derivative of the particle density along the azimuthal coordinate vanishes (∂ φ N = 0), then Eq. 2.1 can be written by substituting: where: At Galactic scale, it is common to assume an azimuthal mean field, B = Bφ, such that, as it follows from Eq. 2.2, D rr = D zz = D ⊥ . In this configuration, the problem reduces to that of isotropic diffusion.
In presence of a vertical component of the mean field, as the case of the GC, with B = Bẑ, we obtain D rr = D ⊥ and D zz = D and, in general, D ⊥ = D .

3D in Cartesian coordinates
In this configuration Eq. 2.1 can be written by substituting: In the present work, we consider only the case in which off-diagonal components of the diffusion tensor, D i =j , are null. Under this condition, the previous equation can be simplified as this: 3 Numerical solution of the transport equation

Discretization over grid
In order to solve the transport equation numerically it is necessary to discretise the equation, i.e. to write it on a discrete grid and transform derivative operators into finite differences. In cylindrical coordinates (2D) we consider a grid with two spatial coordinates (r i , z j ) and one momentum coordinate (p q ); the grid spacing is arbitrary and it may be irregular. In Cartesian coordinates (3D) the spatial grid is instead obtained with three coordinates: CR density at a given position, momentum, and time can be written on this lattice as where n is the time step index. In order to replace the derivatives in the transport equation by their finite difference approximations, we mainly adopt the centred difference scheme for an irregular spaced grid: which gives a truncation error O(∆x 2 ) and O(∆x) for uniform and non-uniform grid respectively.

Iteration scheme
We rewrite 2.1 as a time-dependent equation and we find the stationary solution by evolving an initial condition (IC), N 0 ijk , until it relaxes to an equilibrium solution, N ∞ ijk , for which the time derivative vanishes.
Schematically the transport equation can be now written as: where L is the operator which defines the transport equation. In its discretized version, Eq. 3.1 becomes: where ∆t is the time step, and i is now a unique index over the spatial-energy grid. The algorithm we adopt to evolve the solution of the transport equation 2.1 at each time step is described in Sec. 3.3. The convergence criterion is introduced in Sec. 4.5.

Local One Dimensional (LOD) operator splitting method
A well-known approach to find the solution of a diffusive-advection equation is the operator splitting method.
The basic idea of this algorithm is to consider the transport equation 3.1 as a linear sum of different evolution operators (e.g., radial diffusion, vertical advection, energy loss, ...): and for each of them to find a valid differencing scheme for updating N from timestep n to timestep n + 1, as the operator were the only one on the right-hand side of 3.1. The overall evolution in the time step ∆t is obtained by using separately all the operators in sequence. This specific implementation of the method is known as Local One Dimensional (LOD) operator splitting. Clear advantages of this algorithm are that one can discretise independently the different operators using different methods and different boundary conditions. It also allows to have different time steps for the the different subproblems.
The transport equation in cylindrical coordinates (which is obtained by substituting 2.3 for the CR flux divergence in 2.1), can be conveniently written as the sum of 5 operators:

Crank-Nicolson coefficients
The Crank-Nicolson (CN) method is a convenient discretisation scheme forL i since it is second-order accurate in time and unconditionally stable. According to this scheme, the time derivative is obtained by taking the average of the right-hand side at t and t + ∆t, giving: where n os is the number of active operators. The most general expression ofL can be written aŝ (3.5) and consequently 3.4 becomes: which is a tridiagonal set of simultaneous linear equations that we solve at each timestep to compute N n+1 i once N n i and Q i are given. Tridiagonal systems of linear equations can easily be solved by standard methods like Cholesky decomposition (see Sec. 2.9 in [59]) or LU decomposition (see Sec. 2.4 in [59]). CN coefficients L i , C i , and U i , for the transport operators in the cylindrical symmetric version of the transport equation are reported in Table 1. Boundary conditions (b.c.) are also reported in the same table: We assume N (r = R, |z| = H) = 0, and N to be symmetrical around r = 0. Boundary conditions in momentum are given by: N (p = p max ) = 0 and d dp N p 2 p min = 0 (see also B.4).
The three-dimensional anisotropic version in Cartesian coordinates is easily obtained by discretizing L x and L y similarly to L z with different diffusion coefficients (D xx , D yy , D zz ) and by imposing as boundary conditions: N (|x| = R, |y| = R, |z| = H) = 0. The spatial operator CN coefficients in Cartesian three-dimensional grid are listed in Table 2.
Differently than other operators, we adopt for the energy loss a 2-nd order accurate method which is obtained by integrating in momentum the energy loss term and evaluating the integral with the trapezoidal rule. This approach is described in detail in [45] and implemented in the PICARD code.
Following this idea, we end up with a tridiagonal set of linear equations: that we invert with standard methods to find N n+1 i . In Sec. 4.4, we compare the performances of this scheme with the 1-st order upwind discretisation combined with the CN scheme (see Table 1), as implemented in the previous version of DRAGON.

Validation of the numerical solver
In this section we discuss a complete set of numerical tests aimed at assessing the performances of the evolution equation components. In particular, our intent is to investigate both the convergence and the accuracy of the numerical approach that we use.
In order to test the convergence of the numerical solution, we introduce the concept of residuals. In a steady-state condition, the transport equation can be written as: L( x, p) + Q( x, p) = 0 (4.1) One can then introduce the normalised residual as: When solving the transport equation on a discretized grid, the quantity R corresponds to a matrix over the spatial-energy grid. To identify this matrix with one single value we compute the normalized l2-residual: where N Q is the number of grid points where Q = 0. During the iteration procedure, the residual decreases and finally converges to a constant value which depends on the grid resolution and on the adopted discretisation scheme.
The tests that we discuss in this section are performed under simplified assumptions regarding both the geometry of the Galaxy and the astrophysical quantities. In doing so, we are able to compare the numerical scheme of each transport operator against an analytical solution and this allows us to evaluate the accuracy of the numerical solution.
We quantify such accuracy by means of the relative error rel . This quantity is defined as follows: where N a i is the analytical solution evaluated at the grid point corresponding to the index i, and N m i is the corresponding numerical solution as obtained after convergence has been attained (i.e., after the residual has reached the plateau value).
In the various tests, we evaluate the accuracy of the numerical solution for different values of the time step ∆t and of the grid spacings and we check if the scaling of rel with the grid step is consistent with the scaling expected from the discretisation order.
We derive the analytical solutions for the test cases used in this section in Appendix B.

Spatial diffusion
We first describe the case of a two-dimensional anisotropic and spatially dependent diffusion coefficient with a steady-state source term. The possibility to simulate inhomogeneous diffusion has been a relevant feature of the DRAGON code since its first version. The transport equation adopted for this particular test, together with its analytical solution, is described in B.1. In particular, we study the case of anisotropic diffusion by setting f ≡ D zz /D xx = 0.1.
The left panel of Fig. 1 shows the evolution with time of the squared norm of the residual R 2 , for a fixed value of the timestep ∆t and for different spatial resolutions of the grid. The grid resolution is given in terms of the number of bins along r and z (we take n r = n z = n).
The right panel of Fig. 1 illustrates the impact of the time-step on the residual, for a fixed spatial resolution.  We show in the left panel of Fig. 2 the comparison between the analytical and the numerical solutions at different times, for a fixed spatial and time resolution. In particular, the plot shows the profile along r (and for z = 0) of the two solutions. As it can be seen, the numerical solution reproduces remarkably well the analytical one.
The relative error as a function of the spatial resolution and for different timesteps ∆t is shown in right panel of Fig. 2. The relative error decreases proportionally to the grid step squared, as a result of the discretisation of the operators L z and L r being accurate up to the second order in r and z for a regular grid. We observe that the scaling remains valid down to a minimum resolution, ∆x ∼ √ D xx ∆t, below which the round-off error dominates the truncation one. 6

Advection
We test the advection operator by assuming a Gaussian IC with σ = 500 pc and a constant (both in intensity and direction) advective wind along z with v w = 100 km/s. Differently than the other cases in which we compare with a steady-state solution, this is a typical Initial Value problem where a differential equation is given together with the unknown function in a given point of the solution domain (in this specific case t = 0).
The left panel of Fig. 3 shows the density profile at different times. As expected from the analytical solution (see B.5), at each time the solution corresponds to the rigidly advected IC for ∆z = v w t s , where t s is the simulation time. The maximum relative difference between the IC and the solution obtained at t s = 100 Myr is of ∼ 10 −3 for ∆z = 20 pc and it exhibits a scaling with the grid size as given by the truncation order.  A second test is performed to test specifically the case with discontinuous advective velocity at z = 0 as it is the case of a constant Galactic wind. The analytical solution is given by equation B.25.
We discretize the advective operator with a backward scheme for z > 0 and a forward scheme for z < 0. The CN coefficients at z = 0 are computed by summing the contribution of the z > 0 and z < 0 semi-intervals computed in the forward and backward schemes respectively (see table A).
In the right panel of Fig. 3 we show the relative error of the numerical solution with respect to the analytical one for different grid sizes. The solution with n z = 512 points reproduces the analytical function better than 1%. In the same plot we compare the numerical solution obtained by assuming vanishing CN coefficients at z = 0 as implemented in the previous version of the code.

Momentum diffusion
We consider the reacceleration equation in B.13 with normalization of 5.1 · 10 −16 GeV 2 /s (v A = 50 km/s) and Q 0 = 1. The corresponding analytical solution is obtained in B.4 and the discretization scheme is reported in table A.
Following the same strategy detailed in 4.1, we consider different values of ∆t and different grid sizes; the momentum interval is 0.1 ÷ 10 2 GeV, the number of intervals we consider is 32 ÷ 512. We show in Fig. 4 (left panel) the comparison between the numerical and the analytical solution for ∆t = 1 Myr.
We notice that the time needed to reach convergence is much larger at larger energies (> 10 5 iterations in the case considered in the plot for p > 10 GeV). However, at momenta larger than ∼ 10 GeV the reacceleration operator is usually subdominant when spatial diffusion is also taken into account.  We show in Fig. 4 (right panel) the scaling of the relative error with the number of grid points. We obtain the ∝ n −1 p scaling as expected from the truncation order of the first derivative term in L p .

Energy losses
The analytical solution of the energy-loss equation B.7 is shown in Sec. B.3.
We compare here the two methods presented in Sec. 3.4: a first-order Crank-Nicolson scheme as implemented in the previous version of the code, and a second-order more accurate scheme. Similarly as in the previous sections, we perform the numerical tests considering different values of ∆t and different grid sizes.
We compare in Fig. 5 the analytical solution with the numerical one as obtained for ∆t = 1 Myr with the first-and the second-order discretisation schemes. The reader can appreciate the better accuracy of the second-order approach. With n p = 32 grid points in the given momentum range, the first-order scheme produces a ∼ 25% error with respect to the analytical solution, while the second-order scheme accounts for a ∼ 2.5% error only. We remark that this improvement in accuracy is obtained without a significant increase in the number of iterations required to reach convergence.
Absolute and the relative errors at p = 10 GeV as a function of the grid size is shown in Fig. 6. Expected scalings are correctly reproduced by the numerical solution. Remarkably, the second-order scheme provides a 10 −2 accuracy already with 64 grid points, while the same accuracy is reached by the first order scheme with more than 512 points.

Convergence and variable time step
In a realistic simulation, the code should be able to propagate particles in a huge range of energies and timescales.   As an example, spatial diffusion is determined by a diffusion coefficient whose dependence on the particle momentum is a power-law with slope δ ∼ 0.3 ÷ 0.7 (see Appendix C.8). That implies that high-energy particles can diffuse much faster than the low-energy ones. In particular, if the range of energies under scrutiny spans several orders of magnitude the difference in the diffusion timescales can be quite large. As an example, with δ = 0.5, D 0 ∼ 1.2 × 10 28 cm 2 s −1 and an halo H = 4 kpc, a 10 TeV particle have a diffusion timescale of t d = 2 Myr, which is 100 times faster than a 1 GeV particle. In terms of the numerical solution of the diffusion equation, the code must evolve with a time-step ∆t that has to be sufficiently small to correctly follow the diffusion of high-energy particles and it requires a time to reach convergence which is of the order of the longest timescale at the smallest energy. Because of this, the number of iterations that must be performed in order to find the correct numerical solution through all the energy spectrum could be extremely (and unnecessary) large.
An efficient algorithm to significantly reduce the computational effort has been proposed by [34] and a key feature of GALPROP. The basic idea is to change the value of ∆t within the single run, starting from a large time step and reducing it after a certain number of iterations is performed.
In the example we are considering here, we further assume source distribution having a Gaussian profile along z with characteristic scale 100 pc. We then start our simulation with ∆t max = 100 Myr and we reduce the time step by a factor of 1.1 each m = 5 iterations, until ∆t min = 1 Myr is reached. We remind from Sec. 4.1 that the stability condition for diffusion reads: ∆t/t d < 1, and we choose ∆t as half of the corresponding timescale.
We evaluate the performance of the code under this approach by comparing the local spectrum against the one we obtain by using a constant time-step with ∆t const = ∆t min . The results of the comparison are shown in Fig. 7.
As one can clearly see, by using a constant time-step, the numerical solution reproduces the analytical one only after ∼ 10 3 iterations. In fact, the use of a variable ∆t provides a very efficient way to reduce the number of iterations and an accurate numerical solution is found after 180 iterations only.
By showing the right panel of Fig. 7, we point out how the overall accuracy found by using a variable ∆t is affected by changing the value of m, namely the number of iterations for each time step. Since the total number of iterations is simply proportional to m, a trade-off between accuracy and computational time is required to choose this parameter.
As a second example, we consider the loss equation. We choose a typical value for the energy loss normalisation, b 0 = 10 −16 GeV/s, corresponding to the synchrotron energyloss rate of an electron propagating in a ∼ 5 µG magnetic field. Under these conditions, the associated timescale, t l , ranges from ∼ 3 Gyr at the lowest momentum to ∼ 0.3 Myr at 1 TeV. As in the previous example, we compare the case with a constant ∆t cst = t l,min /2 and a variable ∆ t , starting from the largest timescale and decreasing until ∆t cst is reached. The figure 8, where we show the accuracy with respect to the analytical solution for different choices of the parameter m, confirms a similar behaviour as for the diffusion case.
In conclusion, we found that the variable ∆t scheme provides a large improvement in terms of number of iterations. As pointed out by [45], this procedure could however have a shortcoming. As shown in the above two examples, the numerical scheme parameters, as the largest and the smallest timestep, strongly depends on the physical parameters (e.g., diffusion coefficient, Alfvén speed, etc.). Moreover, the same author discussed how changing the physical parameters in a GALPROP run could generate a final output far from convergence (see also the appendix A.1.1 in [61]). In fact, GALPROP solicits the user to tune the numerical parameters governing the time steps accordingly to the physical model chosen.
On the contrary, DRAGON2 code computes the minimum and maximum time steps according to the active operator timescales. This strategy ensures numerical convergence for most of the physical parameter space.

Physical applications based on new features
In this section we present the capabilities of the new features in DRAGON2 by showing a few example applications. These features are motivated by several pieces of evidence (partially discussed in the following sections) calling for a more realistic description of the CR transport in our Galaxy than the simplified one-zone model implemented in most of the literature.
This chapter is organised as follows.
• In Section 5.1, 5.2, and 5.3 we focus our attention on spatial-dependent and anisotropic diffusion, i.e., the most characterising novelties of DRAGON2. Besides being theoretically motivated, the need for implementing these features comes from γ-ray observations, currently the most effective tracers of the CR distribution in the Milky Way. In particular, we focus here on the so-called gradient problem and the slope problem.
The first anomaly -already noticed in pre-Fermi data (e.g., with COS-B [62]) -consists in a discrepancy between the radial profile of the γ-ray emissivity (inferred by the γray longitudinal profile) and the one computed with conventional propagation models, using an injection term based on pulsar or SNR catalogues: The former appears to be flatter along the Galactic plane at large Galactic radii. This problem was recently confirmed by Fermi-LAT observations [63,64]. Although several possible explanations have been proposed (e.g., a flatter CR source distribution in the outer Galaxy or a strong Galactic wind), this discrepancy may be the signature of a faster perpendicular escape through the halo in the inner region of the Galaxy [27].
The second anomaly is about the γ-ray spectrum: The comprehensive collection of GALPROP-based models provided by the Fermi-LAT collaboration in [65] underpredicts the data above 10 GeV in the inner Galactic plane; this discrepancy is present in all the considered setups, with different prescriptions for the gas and source distributions. This high-energy excess is the signature of a harder CR spectrum in the inner Galactic plane, as shown also in [15,66,67]; according to [68] the hardening can be explained by a progressively harder scaling of the diffusion coefficient in the inner Galaxy (see Appendix C.8).
• In Section 5.4, we investigate the consequences on the proton spectrum of including the hadronic energy losses due to pion production and we quantify the impact on the propagated proton slope.
• In Section 5.5 we briefly discuss the possibility of propagating transient sources. This feature is relevant when CR injection cannot be assumed to be a continuous process in space and time. The imprint of a nearby recent SN event on the locally observed CR spectrum has been recently revised in [69] and [70]. Under peculiar circumstances (in particular, highly anisotropic diffusion along regular magnetic field lines connecting the source to the Solar System), the bulk of low-energy hadronic CRs observed locally could be the result of a single recent event. Motivated by these considerations, we include in our code the possibility to follow the propagation of CRs emitted by time-dependent sources, and present in detail a physical case with anisotropic diffusion.
• As mentioned in Section 3.1, one of the key features of the DRAGON2 code is the possibility to perform CR transport over a non-uniform spatial grid. This property proves useful whenever one wants to study a certain region of the Galaxy (e.g., the Galactic Centre) with much higher resolution than the rest of the domain. In Section 5.6 we present two illustrative examples in which the non-uniform grids can be convenient in terms of computational time.
• The relevance of reacceleration in CR propagation is still unknown: Local CR observables cannot provide strong constraints on this process given the uncertainty on our understanding of ∼ GeV diffusion. Since the new version of the DRAGON2 solver includes an alternative treatment of the boundary conditions in momentum, and a second-order scheme for the energy losses, it is interesting to compare the outcomes with what we would obtain by using the older solver. To this end, we describe the case of leptonic r [kpc ] propagation, where the interplay of diffusion, energy losses (due to synchrotron emission and Inverse Compton), and reacceleration produces a peculiar feature in the spectrum at intermediate energies. We discuss this issue in Section 5.7.

Inhomogeneous diffusion coefficient: The gradient problem
Here we show how to solve the gradient problem in a position-dependent scenario with DRAGON2 working in 2D mode. DRAGON featured already inhomogeneous diffusion, this numerical test is shown to confirm one of the main results obtained with the previous version of the code. Assuming cylindrical symmetry and a purely azimuthal structure of the regular GMF, only the perpendicular diffusion coefficient (D ⊥ ) needs to be considered.
We assume D ⊥ to be spatially correlated to the turbulence strength, hence to the CR source density q(r). Being the exact relation between those quantities poorly known, we parametrise it as in the PropToSourceTerm option described in Appendix C.8, with δ = 0.5 and η = 1. For the source term we choose the Ferriere2001 model, with a slope for momentum of 2.3.
We solve the 2D diffusion equation by adopting a variable time step ∆t, starting from the largest time step ∆t = 16 Myr and reducing it by a factor 0.5 each 100 iterations, until ∆t = 10 −4 Myr is reached.
The main result is reported in Fig. 9: proton radial profile flattens, thus ameliorating the gradient problem, when the value of τ is increased (see also [27]).
We show further tests in Fig. 10, where we focus on the model characterised by τ = 0.7. In left panel we show the effect of different spatial resolutions on the accuracy of the solution. A reference run with a very fine grid (∆r = 60 pc) is used as a comparison. It is interesting to notice that, due to the smooth distributions of all the relevant quantities, a r [kpc] percent-level accuracy all through the Galactic plane can be reached with a grid with ∆r = 1 kpc. In fact, we find that in most of the applications, a spacing of ∼ 200 ÷ 500 pc is adequate to obtain an accurate solution, well below the level of 1%.
In right panel we consider the impact of changing the radial boundary condition. We find that the solution, in the region of interest (r < 10 kpc), is mildly dependent on the position of the boundary, and choosing a value for r larger than 15 kpc provide a negligible improvement in the accuracy of the solution.

Variable scaling of the diffusion coefficient
According to [68], the γ-ray longitudinal profile can be successfully reproduced by a progressively harder scaling of the diffusion coefficient in the inner Galaxy.
In DRAGON2, this scenario can be obtained by using the model VariableSlope for the diffusion coefficient (see Appendix C.8). In [68], a similar result was obtained with DRAGON by adding to the CN coefficients the additional terms containing the spatial derivatives of the rigidity power-law.
The spectra at different locations are shown in Fig. 11. In left panel we reproduce the case as in [68] with a = 0.045 kpc −1 and b = 0.126 (giving δ(r ) = 0.5). These values provide a significant hardening in the inner Galaxy, compatible with the γ-ray data (see [68] for the details) .
The CR proton slope at different radii can be compared with the following prescription: which would be the generalisation of the leaky-box model, and valid in the thin disk approximation.  In right panel of Fig. 11 we compare propagated proton spectra with those expected from Eq. 5.1. In the specific case we are considering, the relation between injected and propagated slopes remains valid at large radii (in particular, it is obtained at the Sun position). As a consequence, a conventional propagation model tuned on local CR observables is still in agreement with local measurements (e.g., B/C) even if the spatially-dependent scaling is introduced.

Anisotropic diffusion coefficient
In this section we show how a radial dependence of the CR spectral index, similar to that discussed in the previous Section, may be obtained in a simplified framework of anisotropic diffusion.
Since CR diffusion is caused by their scattering off the magnetic field perturbations, and since such perturbations are usually expected to have components both in the perpendicular plane and along the parallel direction (with respect to the direction of the regular GMF), CRs can experience both a parallel and a perpendicular diffusion. As discussed in [58], numerical simulations of charged particle propagation in turbulent MFs found that these two components of the diffusion tensor are characterised by a different scalings with respect to the CR momentum.
In order to test such scenario, we consider a GMF with two components: • A purely azimuthal component, lying on the Galactic disk.
• An out-of-plane component, directed along the z-axis and confined within the bulge (R < 2.9 kpc). This model is a simplified version of the Jansson2012 model actually implemented in DRAGON2 and based on [71] (more details in Appendix C.2.1).
We expect therefore diffusion along the r-direction to be purely perpendicular, since the GMF has no radial component, while the diffusion coefficient along the z-axis is given by the sum of a parallel and a perpendicular term: For the normalization of the diffusion coefficients, we assume D 0, = 10 28 cm 2 s −1 and, following [58], D 0, /D 0,⊥ = 30. As already said, parallel and perpendicular diffusion are characterised by a different dependence on particles momentum: following again the results found in [58], we assume that δ = 0.33 and δ ⊥ = 0.5. The profile along r of the diffusion coefficients D r and D z for particles with p = 1 TeV is shown in the left panel of Fig. 12.
Concerning the geometry of the halo, we assume its radius to be 20 kpc, and its height to be 4 kpc. We assume as a source term a Gaussian disk with momentum power-law injection: We show in right panel of Fig. 12 the steady-state energy spectrum on the Galactic plane at different distances from the GC. For low values of R (within the bulge), parallel diffusion dominates and, as a consequence, a significant hardening can be noticed in the propagated  slope; on the other hand, for larger values of R, the slope is steeper, and tends to 2.8, i.e. the value expected for perpendicular diffusion only.

Pion momentum losses
A novelty introduced in DRAGON2 with respect to earlier versions is a new implementation of pion production as a continuous loss term, in addition to ionisation and Coulomb losses.
In our numerical tests we consider a homogeneous source term confined in a disk (with scale height 100 pc); the energy-loss due to pion production relies on the analytical parametrisation reported in Section C.10.6, with the scale height z losses = 100 pc.
Concerning the diffusion coefficient, we refer to the Eq. C.23, and implement a standard diffusive regime corresponding to D 0 = 1.8 · 10 28 cm 2 s −1 at 1 GV, and δ = 0.5; The diffusive halo height is set to H = 4 kpc.
In Fig. 14 we show the results of our numerical tests: We notice that the impact of the pion-production energy loss term on the proton flux ranges from 5% at 100 MeV down to few percent at energies larger than 100 GeV. The impact on the slope is 0.5%.

Anisotropic diffusion from a transient source
The study of a transient source is relevant in many different context (e.g., to describe the Galactic centre activity). For this reason, we show here how DRAGON2 is able to follow the evolution of CRs emitted by an energetic source in a short event. Since the source is pointlike and, in general, far from the centre of the coordinate system, we exploit the 3D mode; moreover, we consider anisotropic propagation, with D r > D z , inspired by a quasi-linear theory scenario dominated by parallel diffusion along the regular magnetic field line on the plane.
Our choice of the parameters is comparable with an average SNR event, as summarised here:  In Fig. 15 we show the evolution of the spectrum on the x axis, at a distance of 1 kpc from the source. For each rigidity the flux peaks after a time t compatible with the order-ofmagnitude estimate given by D(p) = L 2 /(6t). The plot clearly shows the different timescales associated to CR diffusion at different rigidities, with the bulk of high-energy particles ( 100 GeV) arriving at the observer's position after ∼ 0.5 Myr, and low-energy ones ( 100 MeV) after ∼ 10 Myr.

Non-equidistant binning
Sources at the GC (3D mode) We consider a three-dimensional source term given by: where x s = y s = 250 pc and z s = 100 pc represent the size of the source along the x, y and z axes, respectively. We consider a purely diffusive and isotropic scenario characterized by a diffusion coefficient defined as: with D 0 = 10 28 cm 2 s −1 .
We study the diffusion of CRs in this setup by performing a run characterised by a time-step ∆t = 0.01 Myr within a cube of edge 10 kpc. We find that, independently on the spatial grid used in the run, the numerical solution becomes constant after 2×10 4 iterations. We set the spatial resolution along the x and y axes to 156 pc (corresponding to 65 bins in each directions) and we test three different setups for the grid along the z axis: • An Equidistant Binning (EB) with n z = 201, corresponding to a constant resolution of 50 pc from z = −5 kpc to z = 5 kpc.
• A Not Equidistant Binning (NEB) with n z = 27, where the bins width grows from 50 pc for |z| ≤ 100 pc to 100 pc in the interval 100 pc ≤ |z| ≤ 500 pc and then up to 1.36 kpc in the range 500 pc ≤ |z| ≤ 5 kpc (going to larger values of z, the width of each bin is larger than the previous one by a 50%).
• An EB with the same number of bins as the NEB setup that we have just described (27), corresponding to a spatial resolution of 385 pc.
The three different grids in the |z| ≤ 2 kpc region are depicted in left panel of Fig. 16.
The results of the different runs are shown in the right panel of Fig. 16. In particular, the plot shows the CR density profile along the z axis, computed at x = y = 0 pc. As it can be seen, the numerical solution found with NEB is indistinguishable from the one given by EB with n z = 201. A crucial factor that has to be taken into account, however, is that in the former case the solution is found in 9222 seconds, while in the latter the run lasts for 1196 seconds before a stable solution is found. As for EB setup with n z = 27, the runtime is approximately the same as in the NEB case, but the solution found with such a binning can be slightly wrong, in particular in the peak region, where the error can reach 5%.
From these results, one could infer that in a scenario like the one that we considered in our test, which is characterised by a relatively spread CR distribution, using NEB determines a small increase in the accuracy of the numerical solution. Still this improvement comes at no computational cost, making the feature serviceable also in this case.

Sources on a disk (2D mode)
When particles are confined in a region that is much smaller than the size of the Galaxy, NEB can represent a very useful and necessary instrument to obtain a precise solution in a relatively short time. To provide an illustrative example of such a case, we consider a one-dimensional scenario in which we have a Gaussian source term defined as: where z s = 100 pc. As before, we consider a purely diffusive case, but this time we assume that the diffusion coefficient drops by up to three orders of magnitude in the source region: where D 0 = 10 29 cm 2 s −1 . Under a physical point of view, one can motivate this decrease in the diffusion coefficient as a consequence of the stronger turbulence that characterises the region of the source.
We study the propagation of CRs in the -5 kpc ≤ z ≤ 5 kpc region for three different setups of the binning along the z-axis: • An EB with n z = 501, corresponding to a constant resolution of 20 pc.
• A NEB with n z = 31, where the bins width is 20 pc for |z| ≤ 100 pc and then grows to 50 pc, 100 pc, 500 pc and 1 kpc as larger values of z are considered.
• An EB with n z = 31, i.e. the same number of bins of the NEB setup described above. This number of bins corresponds to a spatial resolution of 333 pc.
Left panel of Fig. 17 illustrates the binnings corresponding to the three setups described above in the region |z| ≤ 1 kpc.
The profile along the z-axis of the numerical solutions obtained for the three cases are shown in right panel of Fig.17. As it can be clearly seen, the solution that can be obtained with NEB overlaps perfectly with the one that is found in the EB case with n z = 501, while the solution that characterises the EB with n z = 501 appears to be wrong by more than two orders of magnitude. The advantage of using NEB is here evident, since by going from n z = 501 to n z = 31 the runtime decreases considerably (it goes from 46 seconds to 4 seconds), without any loss in the accuracy of the solution.

Interplay of diffusion, reacceleration, and leptonic energy losses: the role of the boundary condition in momentum.
The lepton spectrum is significantly shaped by the energy losses (∝ E 2 ) due to synchrotron emission and inverse Compton scattering (see C.10 for the relevant formulas). We consider here a realistic setup in which the energy loss term is coupled with other operators: reacceleration and diffusion. A significant reacceleration, in combination with ICS and synchrotron losses, produces a hump-like feature in the spectrum: In this section we aim at characterising this feature in a realistic setup.  As for the hadronic tests, we consider a homogeneous source term confined in a disk (with scale height 200 pc); the energy-loss term is taken as follows: In Fig. 18 we show the CR spectrum at the Sun position, and we compare the new prescriptions for the discretisation of reacceleration and energy loss operators (described in detail in section 3) with the ones implemented in the first version of DRAGON.
The hump is clearly visible in the spectrum in both cases, and give rise to a peak at 2.5 GeV: slightly higher than the energy where the energy-loss and reacceleration timescales are the same (∼ 1 GeV). This feature is caused by the interplay between the two competing effects of energy loss (which is responsible for a downwards flux in momentum space) and reacceleration (which is a diffusive term in momentum space, and is responsible for an upwards flux due to the monotonically decreasing injection spectrum).
The DRAGON solver is characterised by a first-order discretisation for the energy losses, and the boundary conditions in momentum are obtained by imposing null derivatives for the The reference runs are shown in Fig. 18: We set the minimum momentum of CR injection at p inj min = 0.1 GeV and the boundary of the simulation at p min = 0.1 * p inj min . In this setup, the two solvers exhibit a 30% difference at all energies due to the better accuracy of the second-order scheme.
We then compare each solver with the case in which p min = p inj min . We notice that the solution obtained with the new solver is practically unchanged by this prescription, while the solution of the DRAGON solver, with the old boundary condition, is significantly different at lower momenta. We conclude that the momentum boundary condition of DRAGON2 enables us to impose p min = p inj min with negligible effect on the solution accuracy. In Fig. 19, we consider a reference run with n p = 512 points/dex and consider the impact of lowering the resolution. From this plot, we came to the conclusion that a resolution of n p = 16 points/dex (i.e., 64 points in total in the momentum range where the source is injecting particles) yields an error as large as 15% in some portions of the spectrum. Even with a second order discretisation method, a ∼ 1% accuracy in the whole momentum range requires a much larger resolution (n p ∼ 128 points/dex).

Conclusions
The simulation of galactic CR propagation plays an essential role in understanding the properties of the galactic sources and of the interstellar medium. In this paper we introduced the new version of the publicly available galactic CR propagation code DRAGON2.
To interpret the recent measurements of primary CR fluxes and their secondary products, the code was rewritten to solve the complete CR transport equation in 2D and 3D, for both leptonic and hadronic species, under very general assumptions.
The equation solver has been updated for all the operators: anisotropic and inhomogeneous diffusion, advection, reacceleration, and energy losses. As discussed in the text, some of these features have been already exploited with modified version of DRAGON to study diffusion emissions and local spectra (see, e.g., [68,72]), that are now fully tested and included in the public code release. For each operator, we discussed the numerical scheme adopted in DRAGON2 and we tested against the corresponding analytical solution.
In doing so, we provide a complete set of analytical solutions that can be used as a benchmark test for any CR transport framework.
DRAGON2 features a wide collection of up-to-date models for all the relevant astrophysical ingredients involved in the computation: source term, magnetic field, interstellar radiation field, gas; all the distributions implemented in the code are presented in the paper.
We emphasised the novel features of DRAGON2 with respect to its previous version by simulating relevant physical cases where these features play a crucial role. In particular, we showed the implications of a different normalisation and rigidity scaling of the diffusion coefficient in different regions of the Galaxy, we described the case of CR anisotropic propagation from a transient source, we considered the possibility of a refinement of the grid in a peculiar interesting region, we discussed the interplay of diffusion, reacceleration and energy losses and the impact of these processes on the propagated leptonic spectrum.
DRAGON2 is an open source code under GNU general public license 7 distributed as a git repository (see www.dragonproject.org), and it can be easily upgraded or re-used to describe a wider range of physical conditions with respect to those treated in this paper.
Forthcoming publications will focus on additional aspects not covered in this paper, e.g., the general (non-diagonal) treatment of anisotropic diffusion, and a comprehensive description of the cross-section network.

A CN Tables
Operator

B Analytical solutions to the transport equations
In this section, we provide different analytical solutions to the transport equation defined in 2.1. In order to derive these solutions we will consider a series of simplified scenarios in which only one operator of the transport equation plays a role, the others being set to zero. Analytical solutions derived in this way are valuable to test the performances of the numerical solver for each transport operator. As mentioned in Section 4, we can compare the numerical scheme of each operator against its analytical solution, allowing to evaluate the quality of the convergence reached by the corresponding discretisation scheme.

B.1 Two-dimensional inhomogeneous and anisotropic diffusion equation
We provide here an analytical solution for the two-dimensional diffusion equation: in the domain extending from 0 to R in r direction and from −L to L in z direction, and boundary condition N (z = ±L, r = R) = 0 and ∂N/∂r(r = 0) = 0. We assume inhomogeneous and anisotropic diffusion coefficients, namely: where f is the parameter defining the anisotropy between diffusion in the z and in the r direction.
Following [45], we prescribe an analytical solution as: and we find the corresponding source function that fulfils Eq. B.1 as:

B.2 Three-dimensional anisotropic diffusion equation
We consider the equation: which describes a homogeneous and anisotropic three-dimensional diffusion. The diffusion is assumed to be confined in the domain which extends from −L x to L x in x direction, from −L y to L y in y direction and from −L z to L z in z direction. The diffusion coefficients can be written as: where f y and f z are parameters that quantify the anisotropy of the diffusion along the y and z directions. In analogy with the two-dimensional case described above, we prescribe an analytical solution given by: In order to satisfy Eq. B.4, the source term must take the following form:

B.3 Energy-loss equation
Taking only losses and source terms into account, the steady-state equation becomes: We assume different power-laws for both the loss term ratė and for the source term: where p 0 is a reference momentum. The analytic solution of Eq. B.7 can be found by integrating both sides in momentum: which gives: The integration constant C is found by imposing the boundary condition N (p = p max ) = 0:

B.4 Reacceleration equation
In order to find an analytical solution for the reacceleration equation in the momentum range [p min , p max ]: we assume a simplified expression for the momentum diffusion coefficient: while the source term is assumed to be the same as in B.9. The transport equation B.13 becomes: where x ≡ p/p 0 and f ≡ N/p 2 . A first integration of both sides is performed : The integration constant C is given by the boundary condition ∂ ∂x f = 0 for x = x min : A further integration of Eq. B.16 gives: The second integration constant B is derived under the assumption that f = 0 for x = x max : Finally, by substituting the integration constants in Eq. B.18, we obtain:

B.5 Advection equation
The advective transport of N (z, t) by the velocity field v w (z) is described by the equation: The solution of B.21, assuming v w constant along z, is determined by an initial condition N 0 (z): which is just the initial function N 0 shifted by v w t to higher z (for v w > 0) or to smaller z (for v w < 0). In Section 4.2 we assume a Gaussian initial condition centred at z = 0: and we assume positive v w .
In addition, to test the advective transport by the velocity field v(z) = sgn(z)v z , as it is the case of a Galactic wind, we assume a simplified one dimensional equation for N (z): with boundary condition N (|z| = L) = 0. The solution to this equation can be written as: where w can be found by solving Eq. B.24 for z > 0 and it turns out to be: By integrating the Eq. B.24 in the neighbourhood of the galactic disk (i.e., by computing lim →0 − dz . . . ), we end up with:

C Astrophysical ingredients
In this section we will present a review of all those astrophysical ingredients (gas and sources distributions, galactic magnetic fields, ISRF, relevant velocities and energy losses) that play a role in our numerical simulations. These ingredients are not always known with enough accuracy, and different models describing them can be found in literature. It is important to underline that using different models for the same ingredient can lead to systematic uncertainties in the outcome of our numerical simulations and for this reason we implement in the code at least two different models for each ingredient. The modular structure that we introduce in DRAGON2 allow the users to easily extend each module with a new model.
Please note that the names that we associate to every model are the same as those with which these models have been called in the code.

C.1 Gas distributions
Interstellar gas plays an important role in the process of CR energy loss and in secondary production. The radio emission is used to trace its distribution (for a review see [73]). Here we shortly report the different distributions of the interstellar Hydrogen: ionised (HII), atomic (HI), and molecular (H 2 ) phase as implemented in the code. Radial distributions of the two most abundant components (HI and H 2 ) for the models implemented in DRAGON2 are shown in Fig. 20. As customary, we assume that Helium distribution follows the Hydrogen one with a ratio of 0.11 between the number density of the two species.

C.1.1 HI gas distributions
For the neutral Hydrogen we adopt the following models: Gordon1976 This model is based on HI relative distributions as reported by [74] (see their Table 1). The gas density is assumed to be a Gaussian along z with radially binned plane density (n 0,i ), with bin positions r i and half-width at half-maximum (z 1/2,i ). As pointed out by [75], distances inside the solar circle (as r i , z 1/2,i ) scale proportionally to R , while density (as n 0,i ) scales inversely to it. We adjust the model quantities according to the value of R adopted in the simulation, i.e 8.5 kpc.
Nakanishi2003 In [76], a three-dimensional distribution of HI gas in the Galaxy is derived from available HI survey data and the azimuthally averaged surface density and scale height as function of galactocentric radius are shown in their Figures 4 and 10. We fit these quantities with a third-order polynomial (after re-scaling them to R = 8.5 kpc as described in the Appendix D of [75]) and we compute the gas density as: n HI (r, z) = n HI,0 (r) exp − ln 2 z z 1/2 (r) 2 (C.1)

Ferriere2007
We use the analytical fit provided by [77] for r > 3 kpc and based on HI distributions reported in [78] for within the solar circle and in [79] for the galactic region outside it. For the inner Galaxy (r < 3 kpc) we use here an azimuthally averaged version of the inner Galaxy HI three-dimensional distribution given in [80].

Pohl2008
The reconstructed spatial distribution of neutral gas in the Galaxy is obtained by [81], who have kinematically deconvolved the LAB HI survey by [82] with appropriate modifications reflecting the larger single-cloud linewidth and the galactic warp and flaring in the outer Galaxy 8 . In Fig. 21, the gas density pattern along the Galactic Plane (averaged for z) is shown when three-dimensional propagation mode is selected.

C.1.2 H 2 gas distributions
Molecular Hydrogen (H 2 ) is not directly observable and is traced by the 2.6 mm CO emission line. The CO-to-H2 conversion factor is called X CO ; both its normalisation and spatial dependence are affected by large uncertainties.
For modelling the H 2 component we adopt the following distributions (all these models are rescaled for a constant X CO = 10 20 cm −2 K −1 / km s −1 and they are multiplied by the X CO profile chosen among the models described in C. 1

.3):
Bronfman1988 An axisymmetric model for the H 2 is provided by [75], as a fit of the combined Northern and Southern CO surveys performed with the Columbia Millimeter-Wave Telescope. The best-fit model parameters are given in their Tab. 3.

Nakanishi2006
This model is described in [83] and is based on the [84] CO survey. As for the corresponding Nakanishi2003 HI model, we take the radial profiles of the H 2 density at the midplane and of the FWHM from their Tab. 1 and we rescale for R = 8.5 kpc.
Ferriere2007 As for the Ferriere2007 HI model, we implement the distribution as given in [77] for r > 3 kpc and in [80] for r < 3 kpc.

Pohl2008
The same technique utilized for the neutral gas was applied to derive the molecular gas distribution in the Galaxy [81]. The analysis relies on the CO data by [84], while for the inner Galaxy is based on a gas-flow simulation using smoothed-particle hydrodynamics (SPH), introducing a realistic barred gravitational potential derived from the observed COBE/DIRBE near-IR light distribution [85] 9 . In Fig. 21, the gas density pattern along the Galactic Plane (averaged for z) is shown when three-dimensional propagation mode is selected.

C.1.3 X CO models
We then multiply H 2 distribution by the X CO,20 (normalized to 10 20 cm −2 K −1 / km s −1 ) distribution; for the latter we adopt the following radial dependences: Arimoto1996 The authors of [86] studied the radial dependence of the conversion factor by comparing CO luminosities to the virial masses of a large number of giant molecular clouds situated at various galactocentric distances. They obtained the relation (rescaled to r = 8.5 kpc): X CO,20 (r) = 0.9 exp r 7.1 kpc (C.2)

Strong2004
We use the outward radial gradient in X CO adopted in [87] to reconcile the predictions of the GALPROP code, based on the assumption that the distribution of CR sources follows the observed distribution of pulsars, with the γ-ray profiles measured by EGRET/COMPTON.
Ackermann2012 A similar analysis as in [87] was performed by making use of the FERMI data in [65]. The radial-dependent X CO is shown in their Figure 25.
Evoli2012 In [88], a two-zone model for the X CO is used. In this model the X CO inner value, the outer value and the border position can be chosen by the user.

C.1.4 HII gas models
For the ionized Hydrogen component we feature: Cordes1991 In [89] the authors make use of different astronomical measurements ( pulsars emission, dispersion measurements and radio wave scattering measurement) to derive a cylindrically symmetric model for the Galactic distribution of free electrons. From this distribution is possible to infer the properties of the ISM ionized component. NE2001 The NE2001 model is described in [90,91] 10 . The density of the Galactic free electrons has been derived from the combined fit of pulsar dispersion measures, temporal and angular broadening of radio pulses, and emission measures. The NE2001 model contains several components, and we consider the most relevant on a Galactic scale: (a) a thick disk with large Galactocentric scale height; (b) a thin (∼ 140 pc), annular disk in the inner Galaxy; (c) a Galactic centre component.
We additionally consider a correction proposed in [92]. By excluding sight lines at low Galactic latitude that are contaminated by HII regions and spiral arms, they find that the height of the thick disk roughly doubles to ∼ 2 kpc.
Ferriere2007 As for the Ferriere2007 HI model, we implement the distribution as given in [77] for r > 3 kpc and in [80] for r < 3 kpc.

C.2 Magnetic fields models
Galactic magnetic field plays a crucial role in lepton energy losses (see Section C.10.1), especially by synchrotron emission. Moreover, as discussed in the Introduction, it enters to determine the diffusive properties of charged CR.
There are several methods to measure and constrain the intensity and the orientation of this component: Zeeman splitting observations [93], infrared, synchrotron and starlight polarisation studies [94][95][96], and measurements of the Faraday rotation measures of Galactic and extragalactic sources [97,98].
The Galactic magnetic field B is usually described as a sum of two components: a largescale regular B reg and a small-scale random part B ran both having a strength of the order of µG in the Solar system neighbourhood (e.g., [99]). Away from this local region, however, the properties of those components are poorly known and one has to invoke phenomenological models.
We consider models for the regular and random component separately.

C.2.1 Regular component
Regular magnetic field is divided in a disk, B disk extending few hundred parsec away from the GP, a thicker halo, B halo , and a less intense vertical -poloidal or X-shaped -components. We implement in DRAGON2 five different models for the regular Galactic magnetic field. All those models are parametrised in cylindrical coordinates (r, θ, z): There are two classes of models, depending on whether the direction of the field in two adjacent arms is the same (axisymmetric, or ASS model) or opposite (bisymmetric, or BSS model).
Sun2007 ASS/BSS In this model the disk field takes the form of a coplanar and constant pitch angle spiral field with components [100] where p is the pitch angle and D 1 (r, θ, z) parametrizes its spatial variation: decreasing exponential behaviour are assumed both for increasing r and z. Its general expression is reported in Eq. (7) of [100]. D 2 (r, θ, z) includes possible reversals and asymmetries in the disk magnetic field distribution. We consider the ASS and BSS options corresponding to the expressions in Eq.s (8) and (9) of [100] respectively. The expression for the halo field, corresponding to a double-torus configuration with a scale-height of few kpc's, is given in Eq. (10) of the same paper.
No vertical component is present in these models.
Pshirkov2011 ASS/BSS The model adopts the same parametrisation as the Sun2007 model. The best-fit parameters, derived from a wider set of Faraday Rotation measurements, is however different as reported in Tab. 3 of [98] for both the ASS and BSS configurations.
Jansson2012 In [71], the authors use a Galactocentric (r, φ, z), and a right-handed Cartesian (x, y, z) coordinate system, with the Sun along the negative x-axis, at x = −8.5 kpc. The Galactic north is in the positive z-axis, and φ = 0 is along the positive x-axis. The magnetic field is set to zero for r > 20 kpc, and in a 1 kpc radius sphere centred on the Galactic centre 11 .
The model features a disk, a toroidal, and an out-of-plane component. 11 The public C++ code has been excerpted from: http://www.ccppastroparticle.com/projects/jf12/. • The disk component is a generalisation of the model derived in [101].
• The toroidal component has a purely toroidal, i.e., azimuthal, component with different radial and vertical extension in the northern and southern regions.
• The out-of-plane component is poloidal and axisymmetric, and specified -at any position (r, z) -in terms of r p (i.e. the radius at which the field line passing through (r, z) crosses the mid-plane (z = 0)). Motivated by the X-shaped field structures seen in radio observations of external, edge-on galaxies, the authors refer to this field as the "X-field" component (see also Fig. 22).
The model also features a striated random (or ordered random) field which is aligned along some particular axis over a larger scale, but whose strength and sign varies chaotically on a small scale. In [71] the energy density of this component is assumed to be proportional to that of the regular field. We refer to that paper for all the details. The relevant plots are shown in Fig. 22.

C.2.2 Random component
The random component of the Galactic magnetic field plays a crucial role in CR diffusion and is expected to determine the diffusion tensor and the size of the diffusion region. Unfortunately, very little is known so far about its properties and its relation with the regular field. From the observational point-of-view, one of the most relevant information comes from the spatial correlation of Faraday Rotation measurements. Using extra-Galactic sources the authors of [103] found a structure function compatible with a Kolmogorov power spectrum up to a scale of ∼few pc with a rms amplitude of ∼ 1 µG. On larger scales a flatter spectrum was found on the basis of pulsar rotation and dispersion measurements with an amplitude 6.1 ± 0.5 µG in the Galactic plane region [104]. Interestingly this value is in good agreement with that inferred from the intensity of the synchrotron emission of the Galaxy (which was found to be dominated by the effect of random field) [105,106].
Synchrotron emission can also be used to estimate the vertical extent (along z coordinate) of the turbulent field, hence that of the cosmic ray diffuse halo. Recent studies [31,107], showed that this quantity must be larger than 2 kpc favouring thicker (L ∼ 10 kpc) halos. The shape of the vertical profile, however, is poorly constrained: exponential, Gaussian or step-like profiles are all compatible with presently available data.
Although disfavoured from those results a correlation between the turbulent and regular fields can not, in principle, be excluded.
For this reasons, in DRAGON2 we implement for the turbulent field the following three different prescriptions: ProportionalToRegular In this model, B rms is given by an user-provided rescaling factor of the regular magnetic field intensity | B reg |.
ExponentialModel This amounts to a double exponential in r and z: This structure resembles that of the simplest (azimuthally symmetric) regular magnetic fields models (see, e.g., [97]).
It was shown that for B 0 = 6.1 µG, R B 6 kpc and z B ∼ 2 kpc this model provides a reasonable description of the synchrotron emission of the Galaxy [105].

C.3 Interstellar radiation field
The InterStellar Radiation Field (ISRF) is relevant for the leptonic energy losses via Inverse Compton scattering (ICS). The low-energy photons involved in this process originate from stars, and are further reprocessed by Galactic dust; CMB photons also contribute with a comparable energy density. We implement in DRAGON2 the following three models: ν [Hz]

Porter2006
With this model, we make use of the public dataset provided with GALPROP [65,111] and available for download on the code website 12 . This file, distributed in FITS format, contains the ISRF spatial distribution for 56 log-spaced frequencies in the range ∼ 0.1 ÷ 36000 µm.
This distribution is based on a realistic modelling of the interactions between starlight and the interstellar matter, which takes into account an accurate knowledge of both the star, gas and dust in the Galaxy, and the infrared emissivities (and spectral shape) per dust grain: a detailed calculation is presented in [105] and updated in [112].

Delahaye2010
We implement the ISRF model described in [113]. In this model, the spectral shape is obtained as a superposition of 6 Planck distributions, with effective temperatures for the CMB, stellar, and ultraviolet components as in table 3. However, these parameters can be also modified by the user.

Vernetto2016
In this more recent result [114], a precise description of the Galactic dust emission is computed following the parametrisation suggested by [115]. The energy density of the radiation fields in the solar neighbourhood is found to reproduce well the measurements of COBE-FIRAS and the sky maps of COBE-DIRBE in the wavelength range 60 ÷ 600 µm. The spectral density as a function of the Galactic position in cylindrical coordinates is provided by the authors as a FITS file.
We plot in Fig. 23 Table 3. Temperatures and relative normalisations for the Planck distributions as implemented in the Delahaye2010 ISRF (see also Fig. 23).

C.4 Source term
Assuming SNR as the main CR sources and with universal injection spectrum, the source term can be written as it follows: where Q 0 is a normalisation factor, q SN (r, z) is the rate per unit of volume, Φ(ρ) is the injection spectrum and ρ c the rigidity cut-off. Moreover, the code allows the user to provid a generic function of space and momentum to be implemented as a source term.

C.4.1 Source profile
For the axisymmetric CR source distribution, we adopt the parametrisation suggested by the authors in [116] to model the radial distribution of the pulsar progenitors in the Milky Way, and described by the following function: where A, a, b, R 1 and z 0 are model parameters. In particular, A is a normalisation constant and it is computed to account for the SN galactic rate. R 1 is included to obtain a nonzero surface density at r = 0, which may be not realistic (see also [117]). For what concerns the vertical distribution, the z−dependence in Eq. C.9 reflects instead the assumption of confinement of the sources in the disc. In Tab. 4 we summarize the best-fit parameters characterising the source models based on Eq. C.9, and included in DRAGON2 code as possible source term (see also Fig. 24).
A different parametrisation for the SNR spatial profile is given in [73]. This model traces both type-II and type-I SNR distributions. In particular, the spatial distribution of type-II SNe is traced by the H II regions -produced by their OB progenitor stars -or PSR, instead the type-I SNe follow the distribution of old disc stars.
The SNR profile in this Ferriere2001 model is then given by the sum of two contributions:

C.4.2 Source spectra
The injection spectrum can be chosen between two different models: • a multiple broken power-law with an exponential cutoff: where Θ(x) is the Heaviside function and In doing so, α i , ρ b,i , ρ c,i , f i are free parameters provided by the user and can be different for each species i.
• a log parabolic power-law: with free parameters:

C.4.3 Source normalization
We foresee that, in some applications, the source term requires to be normalised. The user is allowed to provide the galactic source luminosity, L SN , such that the source term is normalised as it follows:

C.5 Wind velocity
The wind velocity is important for adiabatic energy losses and advective transport.
In the new DRAGON2 code we adopt three different parameterisations:  Figure 25. The spiral patterns Wainscoat1992 [123] and Steiman2010 [124]. The Solar System location is identified with a yellow circle.

C.6 Alfvén velocity
For Alfvén velocity, entering in momentum diffusion calculations, the user can choose to provide a spatially constant value (Constant) or to consistently calculate v A as function of position in terms of the magnetic field and ionised gas density distribution (SpatialDependent): In the latter case, we assume that v A is bounded by 100 km/s even for very small gas densities.

C.7 Spiral galactic pattern
We introduce in the code the possibility to superimpose a spatial pattern to the distributions of different astrophysical quantities: source term, gas, ISRF and magnetic field.
In order to apply a generic pattern f ( x) to the distribution c( x), we first compute the new distribution: Then, we rescale the altered term c imposing that the volume integral is conserved Several astronomical observations point to the presence of a spiral arm structure in our Galaxy. There has been a general consensus on the existence of this pattern between ∼ 3 and ∼ 10 kpc, although some disagreement remains on the details regarding the number of spiral arms and their geometry.
Many astronomical data may be used to trace the details of the spiral structure. The arms are generally modelled as logarithmic spirals, and two-, three-, four-arms patterns can be also found in literature.  Table 5. Parameters for the spiral pattern models in Fig. 25. Both models depict four-arms structures. The spiral width is assumed to be 0.2 kpc.

Model name
In particular, the spiral structure may be derived: • From the catalogues of point sources: They reveal the presence of the pattern and permit to fit the arm parameters. For example, the IRAS catalogue of sources in the the 8-25 micron band (mainly red giant stars) was used in [123] to derive a physically realistic representations of the Galactic disk, bulge, stellar halo, spiral arms. The arm parameters derived in this work are adopted in several more recent works [117,125].
• From the study of the far-infrared cooling lines of the interstellar gas. Since the spiral arms, compared to inter-arm regions, show higher gas column densities and star formation rates, they have enhanced line radiation whose observations can be used to trace the arms. For example in [124]  In DRAGON2 we implement the spiral structure with the following formalism. Each arm is modelled as a logarithmic spiral with arm width σ a . The locus corresponding to this curve -in polar coordinates -is defined by: The parameter K is related to the pitch angle p though the following equation: p = arctan(1/K).
The user is allowed to implement an arbitrary number of spiral arms, characterised by the values of K, R 0 , θ 0 , and σ a . In Fig. 25 we provide two examples of spiral patterns, with σ a = 200 pc and arm parameters reported in Tab. 5.

C.8 Spatial diffusion coefficient
As mentioned in Section 2, the interaction between charged CRs and the Alfvén waves travelling through the turbulent Galactic plasma is responsible for the CR random walk, which is well described in terms of a position-dependent diffusion tensor.
We remark that accounting for a spatially inhomogeneous diffusion coefficient is a innovative feature of DRAGON since its first release. Moreover, as detailed below, this feature is essential to interpret recent measurements of CR fluxes and gamma-ray diffuse emissions.
OneZoneModel A constant diffusion coefficient can be adopted to reproduce local observables. It can be parametrised as: where D 0 is the normalisation of the diffusion coefficient, ρ 0 is a reference rigidity, δ is the slope and η allows us to take into account physics taking place at low energy, including some non-linear phenomena such as the dissipation of magneto-hydrodynamics waves by their resonant interaction with CRs [126].

ExponentialHalo
The assumption of diffusion inversely proportional to the magnetic field turbulent component turns into a decrease of D in the halo.
To account for this, we model D as: where z t is a caracteristic scale, and the halo size is ∼ 3z t .
PropToSourceTerm This model was introduced in [88] to provide a solution to the longstanding gradient problem in the outer Galaxy. The diffusion coefficient follows the radial source profile as: where τ is a free parameter.
VariableSlope This model was introduced in [68] in order to reproduce γ-ray emissivity in various Galactic regions. The slope in Eq. C.23 is replaced by a function of the Galactic radius: with local normalization δ(r ) = 0.5. A is the free parameter which can be fitted against γ-ray data. Authors in [68] provide as best fit parameter A = 0.035 kpc −1 .

C.9 Momentum diffusion coefficient
The turbulent fields affect particles trajectory essentially in two ways: particles experience parallel and perpendicular scattering with respect to the background field (mainly due to the turbulent magnetic fields), but they also experience momentum diffusion or stochastic acceleration. In particular, the nonlinear interplay between particles and turbulent waves/modes is a stochastic process that drains energy from plasma turbulence to particles. In this section, we will focus on stochastic reacceleration of CRs due to resonant interaction with turbulence, in particular with low-frequency waves: we refer to two different kind of waves, namely Alfvénic and Magnetosonic waves.
Alfvénic With this model we adopt the momentum diffusion coefficient for the ISM described in [51]. In doing so: where w is the turbulence level (usually assumed to be ∼ 1), v A is the Alfvén speed and D is the direction averaged spatial diffusion coefficient at each position.
Magnetosonic The diffusion coefficient in momentum for scattering by fast magnetosonic waves is [127]: 28) where L is the scale of turbulence injection, v F is the magnetosonic phase velocity, 1/k d is the dissipation scale, and W (k) = W 0 L −2/3 k −11/3 is the energy density power spectrum.
This model could be relevant in highly-turbulent environments, as exploited in [128] to study stochastic reacceleration inside the Fermi bubbles.

C.10 Energy losses
Nucleons lose energy via ionisation, Coulomb interactions with the interstellar gas and pion production, while electrons and positrons lose energy via inverse Compton, synchrotron emission and bremsstrahlung as well.
In the following, we review how these processes are implemented in DRAGON2. Figure 27 present a comparison of the energy loss time for these mechanisms.

C.10.1 Synchrotron radiation
Synchrotron radiation is emitted by an ultra-relativistic particle interacting with a magnetic field. The emitted power is ∝ m −4 , therefore this process is relevant only for leptons, due to their relatively small mass. A detailed description of the synchrotron radiation theory can be found in [129][130][131].
We implement in DRAGON2 the average loss rate for relativistic electrons in an isotropic magnetic field, as derived in [132] and [55]: where σ T is the Thomson cross section and U B = | B r | 2 + | B t | 2 /8π is the total energy density of the magnetic field.

C.10.2 Inverse Compton Scattering
Inverse Compton scattering (ICS) takes place when highly relativistic electrons and positrons scatter with low energy photon fields. This process can be neglected when treating hadronic CRs.
Because of this process leptons lose a part of their energy and photons are up-scattered into high-energetic γ rays.
We take into account highly energetic e ± CRs which scatter against the "sea" of photons constituted by the Galactic ISRF. The energy-loss rate for an electron with an initial energy i and a final energy f (after ICS) is expressed in the following way [131,133]: The term dN coll dtd i d f is the collision rate, given by where the normalizations N a and temperatures T a of the ISRF components for the local position in the Galaxy are reported in Tab. 3. Fig. 26 shows also a comparison between synchrotron and ICS energy losses: the contributions to ICS losses due to different components of ISRF are shown.

C.10.3 Bremsstrahlung
This process, also called braking radiation, occurs when an electron or a positron is accelerated by the electric fields associated with interstellar ions or nuclei. Bremsstrahlung due to hadronic CRs is then neglected, since because of the heavy masses involved accelerations are very small.
In DRAGON2 we implement the following expressions. For a complete derivation, see [132] and [55].
• Ionised gas and weakly-shielded (WS) neutral gas. This expression applies to the interaction with ionised gas at all energies. In the case of interaction with neutral gas, the formula only applies for γ < 100, also denominated weak-shielded regime on neutral gas.
n He cm −3 (C.35) where M i is the atomic mass, and T i is the radiation length with T H ≈ 62.8 g/cm 2 for H and T He ≈ 93.1 g/cm 2 for He. This expression holds for neutral gas in strong shielding case (γ ≥ 800).
• Intermediate-shielded neutral gas. This regime holds for 100 < γ < 800 and it's given by a linear interpolation of the strong and weak shielding in Eqs. C.35 and C.34 [132].

C.10.4 Ionisation Losses
Relativistic charged particles moving through a material medium interact with electrons belonging to atoms in that same material: the interaction thus excites or ionises the atoms.
This process applies to both hadrons and leptons.
• For ionisation losses suffered by hadrons, a general formula derived from [134] was implemented in DRAGON2: Here,Ĩ s denotes the geometric mean of all ionisation and excitation potentials of the atom: in DRAGON2, we usedĨ H = 19 eV andĨ He = 44 eV, as derived in [134]. Also, q max is the largest possible energy transfer from the incident particle to the atomic electron, as defined by kinematics: q max ∼ 2m e c 2 β 2 γ 2 1 + [2γm e /M ] (C. 38) where M m e is the nucleon mass.
• If the incoming relativistic particle is an electron Eq. C.36 is no longer valid and the ionisation losses is given by [133]: where I i is the ionisation potential, for which we use I H = 13.6 eV, I He = 24.59 eV, and E max = γ 2 m 2 c 2 1+γ is the maximum kinetic energy which can be transferred to the stationary electron. A simplified expression of Eq. C.39 is: Eq. C.40 is different with respect to what is included in GALPROP, where the Bethe-Bloch formula as taken from [132] is considered. However, the differences between the two implementations are always smaller than 1%.

C.10.5 Coulomb scattering
Coulomb collisions in a completely ionised plasma are dominated by scattering off the thermal electrons.
• The corresponding expression for hadronic particles is given by [134]: Here, T e is the electron temperature, and the functional form of W e is given by: The Coulomb logarithm, ln Λ in the cold plasma limit has been derived, e.g., in [135]:  CR hadrons can lose energy through the production of pions that might follow a collision against a nucleus of the interstellar medium. This energy loss has been described for the first time in [134] and it has recently been revisited in [136], where the impact of the new parameterisations for the pion production cross section derived in [137,138] has been discussed.
In our implementation of the energy loss due to pion production, we use the analytical formula provided by [136] which reads: where E is the energy of the proton, while n gas = n HI + 2n H 2 denotes the interstellar gas density. Analogously to [136], to model the energy loss by heavier nuclei we assume that the loss rate increases by a factor A 0. 79 . We notice that the formula above can be applied only for protons with E 1 GeV. At lower energies however this loss mechanism is usually subdominant with respect to the other terms. Helium number density, assumed to be 0.11n H p particle momentum per nucleon Φ CR injection spectrum by SNR q SN rate of CR injected by SNR per unit of volume r radial coordinate in cylindrical coordinate system r e classical electron radius R Sun position with respect to the galaxy centre R galaxy radius ρ particle rigidity (≡ p/Z) per nucleon σ T Thomson cross section U B magnetic energy density v A Alfvén velocity v w Galactic wind velocity w turbulence level x coordinate in Cartesian system y coordinate in Cartesian system z vertical coordinate in Cartesian system or in cylindrical system z t vertical scale-length of diffusion coefficient Z atomic number