Regimes of cosmic-ray diffusion in Galactic turbulence

Cosmic-ray transport in astrophysical environments is often dominated by the diffusion of particles in a magnetic field composed of both a turbulent and a mean component. This process, which is two-fold turbulent mixing in that the particle motion is stochastic with respect to the field lines, needs to be understood in order to properly model cosmic-ray signatures. One of the most important aspects in the modeling of cosmic-ray diffusion is that fully resonant scattering, the most effective such process, is only possible if the wave spectrum covers the entire range of propagation angles. By taking the wave spectrum boundaries into account, we quantify cosmic-ray diffusion parallel and perpendicular to the guide field direction at turbulence levels above 5% of the total magnetic field. We apply our results of the parallel and perpendicular diffusion coefficient to the Milky Way. We show that simple purely diffusive transport is in conflict with observations of the inner Galaxy, but that just by taking a Galactic wind into account, data can be matched in the central 5 kpc zone. Further comparison shows that the outer Galaxy at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$>5$$\end{document}>5 kpc, on the other hand, should be dominated by perpendicular diffusion, likely changing to parallel diffusion at the outermost radii of the Milky Way.


Introduction
The origin of cosmic rays has been the subject of much research since the first detection of cosmic rays in 1912, see e.g. [1] for a review. Baade & Zwicky's [2] hypothesis that supernovae are the most likely energy source of cosmic rays has been strengthened, and significant advances have been made in understanding particle acceleration in supernova remnant environments [1]. Equally important is a thorough understanding of cosmic-ray transport and interactions. In this paper, the focus lies on the understanding of the cosmic-ray diffusion process, which determines the evolution of the density of cosmic rays in a spatially limited environment like the Milky Way. In general, turbulent mixing of collisionless charged particles in magnetized plasma offers interesting challenges not encountered in hydrodynamic turbulent mixing, and thus is highly appropriate to this topical collection. Due to the high electrical conductivity of plasma, the magnetic field itself is mixed by turbulence, and the field lines may become stochastic [3]. To the extent that cosmic rays, being electrically charged particles, follow the field lines, they are spatially mixed the same way the field lines are. In addition, because charged particle orbits are helical, with an energy-dependent size, and because cosmic rays are scattered by orbit scale magnetic fluctuations, they both cross field lines and reverse the direction of motion along them. Moreover, the problem is nonlinear in the sense that the cosmic rays modify their environment by heating, pressurizing, and transferring momentum to the thermal background gas, and by feeding back on the rate of scattering itself. Because the time scales for all these processes are similar to the time scales on which the entire astrophysical system evolves, the problem is inherently nonequilibrium. Developing a full description of the resulting spatial mixing process is one of the most challenging problems in plasma physics, and is necessary not only to understand the transport of cosmic rays, but also particle confinement in laboratory plasma devices, and thermal conductivity and viscosity in both lab and natural plasmas.
Diffusion is usually incorporated through the Parker transport equation Here, is an advection speed, p is the momentum, n is the particle distribution, ̂ is the spatial diffusion tensor, and pp is the momentum diffusion scalar. In this equation, we have assumed isotropy in momentum space. While in the Parker transport equation, sources and sinks of cosmic rays are included within the term S, all losses due to interactions, such as synchrotron radiation or spallation of nuclei on the ambient interstellar medium, are neglected.
Cosmic-ray diffusion is believed to be the dominant process for the transport of cosmic rays in many astrophysical environments [4][5][6][7]. However, the components of the spatial diffusion tensor are also the most elusive parameters in Eq. (1), because they depend on a variety of physical processes: Depending on the ratio of the gyroradius to the correlation length l c of the turbulence, certain effects dominate in the diffusion process, leading to five different transport regimes, see [12]. In Fig. 1, this is shown schematically based on the example of a charged particle that moves through different magnetic field configurations. Along the trajectory, the magnetic field strength between the boxes changes, so that four transport regimes (the 5th regime is a transition between two other regimes) are illustrated. The diagram illustrates how the motion of the particle in the quasi-ballistic regime (QBR) is largely unaffected by the specifics of the field lines as a result of the large gyroradius. The decrease of the gyroradius results in a stronger influence of the magnetic field lines on the diffusion properties of the charged particles, because they are tied closely to the lines, but also because the effects of the cross-field motion become more important up to a certain level. The diffusivity of the field line itself adds to the diffusivity of the following particles due to random walk. Overall, resonant scattering dominates the diffusion process, especially in the RSR. Resonant scattering of particles with a pitch angle Θ 0 at fluctuations of wavelength l is present only when the resonant criterion is fulfilled where = cos Θ 0 is defined as the cosine of the pitch angle. To what extent mirror effects towards lower energies compensate the absence of resonant scattering and lead to diffusion remains a subject of active research [8][9][10][11]. Within a magnetic field tot , which is composed of a coherent component along the x-direction = B x and a turbulent component , the diffusion tensor can be written in matrix form [1,[16][17][18] Transport along the coherent magnetic field is described by the parallel diffusion coefficient ∥ ≡ xx and is complemented by the perpendicular ⟂ ≡ yy = zz . 1 The antisymmetric diffusion coefficients A are often negligible or absorbed into potential drift terms [19]. The exact form of the diffusion tensor ̂ strongly depends on the exact realization of the turbulent magnetic field, including its power spectrum G(k) and the turbulence ratio = b∕B . Some of the parameters can be deduced by comparing transport model predictions to observational data, such as the boron-to-carbon ratio [20], by numerical test particle simulations as presented here, and theoretical models [21].
In particular, the so-called leaky box model [22] of the Milky Way predicts that the cosmic-ray energy spectrum observed at Earth is steepened during propagation: Using the transport equation (1) and reducing it to a stationary case ( n∕ t ≈ 0 ) in which diffusion is the dominant process, we obtain a simplified equation Here, we approximate the diffusion process via the time 1 In general, the two perpendicular directions can be different, yy ≠ zz . However, within the Galactic magnetic field, they are usually degenerate, yy = zz = ⟂ . It should also be noted that perpendicular diffusion does not necessarily describe crossing magnetic field lines.
where d is the escape distance and is the diffusion scalar. Thus, the cosmic-ray density is given by the ratio of the source spectrum and the diffusion coefficient, Assuming diffusive shock acceleration, one arrives at a power-law spectrum at injection, i.e. S(E) ∝ E − s . The diffusion coefficient in quasi-linear theory (QLT) with a wave spectrum that follows a power law also becomes a power law (see e.g. [18]), i (E) ∝ E i with i = ∥, ⟂ . Thus, in case of diffusion-dominance, the energy spectrum of cosmic rays after propagation is steepened, i.e. n(E) ∝ E − s − i [23]. These arguments are based on QLT, in which only linear terms in the distortions of the electromagnetic fields and particle population with respect to the undisturbed fields are being considered to simplify the kinetic equations and to facilitate their analytical treatment. In particular, if we assume isotropic turbulence, we obtain an isotropic wavevector spectrum of the form G(k) ∝ E − , with = 5∕3 [24] for Kolmogorov-like turbulence and = 3∕2 for the Kraichnan type [25]. In QLT, this leads to a parallel diffusion coefficient ∥ ∝ E 2− [26]. While both spectral indices describe the observed solar wind turbulence inertial range within the uncertainties [27], we only consider Kolmogorov turbulence in the following. Based on this QLT prediction, there is an expected difference of about 0.2 in the index of the energy dependence of the diffusion coefficients between the two turbulence models, which results in an underlying uncertainty in our results due to the uncertainties in the turbulence model. Fig. 1 Schematic picture of a charged particle traveling through different magnetic field strengths. The magnetic field strength increases starting from the upper left panel and proceeding clockwise, resulting in an accompanying decrease of the particle gyroradius. The ratio of gyroradius and correlation length determines the predominant diffusion regime of the particle. The possible regimes are the quasi-ballistic regime (QBR), the resonant-scattering regime (RSR), the mirror regime (MR), and the non-resonant-scatter regime (NRSR) [12]. Equivalently, the different regimes can also be represented schematically by changing particle energy between the boxes at constant magnetic field strengths. In each regime, different processes dominate diffusion, such as FLRW [13] and the transport across magnetic field lines [14]. The condition for diffusion in QBR follows from [15] and is generalized here for the inclusion of background fields. n stands for the number of boxes with side length l c . The box and thus the scale that would be necessary to allow diffusion exceeds the width of the region shown in the panel, as nl c ≳ r g ∕l c ≫ r g . The mean-free path is approximately given by nl c , for the smallest n which satisfies the condition for diffusion Previous numerical simulations of particle transport in a combined turbulent field b and homogeneous field B were in agreement with this spectral index up to high b∕B ∼ 1 [16,28,29]. Such simulations are typically performed on a 3-dimensional Cartesian grid with spacing s and number of grid points N grid , injecting relativistic particles from a point source in a homogeneous ( plus turbulent magnetic field in order to derive the diffusion coefficient. Radiative losses are not considered, since they are not relevant here. Recent studies have, however, pointed out that these results need to be interpreted carefully, as the energy range in which the simulations are fully diffusive is strongly limited [10,12,17,[30][31][32]. Using simulations in the fully resonant scattering regime, we quantitatively investigate the energy behavior of the diffusion coefficient as a function of the turbulence ratio b/B. Finally, we interpret these results in the context of recent measurements of the diffuse cosmic-ray flux in the Milky Way.

Simulations of the diffusion coefficients
Our simulations are based on the Taylor-Green-Kubo (TGK) formalism [33], see e.g. [15-17, 28-31, 34]. Building on work that focused on specific parameter points and resolutions, we conduct a systematic simulation and analysis setup, similar to [12] that reveals key conditions on numerical simulation requirements. This method uses the fact that the fundamental solution of the diffusion problem, ii Δf (x i , t) = f (x i , t)∕ t , is known to be a Gaussian whose width is described by the diffusion coefficient ii . The second moment of the deviation in x i provides an analytical solution ⟨ Δx 2 i ⟩ = 2 t ii . The left-hand side of this equation can be calculated in simulations of particles that are emitted from a point source, which is placed in a field composed of a homogeneous component = B x and a turbulent component . Here, we use a Kolmogorov-type spectrum, i.e., isotropic and without intermittency with = 5∕3 , k min = 2 ∕l max , and k max = 2 ∕l min , where l min is defined as the smallest numerically resolved wavelength, and l max represents the largest wavelength used in the simulation. Evaluation of the correlation length in the limit l min ∕l max ≪ 1 yields l c ≈ l max ∕5 [35]. The synthetic isotropic three-dimensional turbulent magnetic field is generated and stored discretely on a regular, three-dimensional Cartesian grid with N 3 grid grid points and isotropic spacing s spacing using the inverse Fourier transform of field vectors ( ) that are computed on a regular grid in threedimensional wavenumber space. Linear interpolation yields the magnetic field at an arbitrary trajectory position between grid points. For discrete step sizes s step = v Δt , the diffusion coefficient can be calculated as i indicates the three spatial directions. Here, we consider parallel diffusion ∥ , as well as perpendicular diffusion ⟂ . Furthermore, t = n Δt is the time after n time steps. We propagate particles on a grid with N grid = 1024 , s = 0.85 pc, l min = 1.7 pc, l max = 82.45 pc, and a step size s step = min(r g ∕5, l max ∕20) to ensure that the gyration motions as well as the fluctuations are sufficiently resolved. With the large grid and correspondingly broad spectrum, particles always find waves for interactions [36]. The comparatively large l grid = N grid ⋅ s ≈ 50l c is chosen in order to reduce problems caused by the continuation of the turbulence at the grid-box boundaries.
While we adopt inertial range Kolmogorov scaling for the turbulence, this is an idealization. Interstellar turbulence is driven on many scales, from as much as 100 pc by superbubbles to kinetic scales by cosmic rays themselves [37]. Moreover, Alfvénic turbulence, unlike hydrodynamic Kolmogorov turbulence is known to be highly anisotropic. While compressibility effects can generate an isotropic component [37], this mechanism is unlikely to result in a simple inertial cascade. Whereas the turbulencedependent energy scaling, which will be determined in the next section, can be scaled from the large scales considered here to small scales straightforwardly, the drivers of turbulence relevant on every scale considered must be taken into account for the normalization of the diffusion coefficients.
For Gaussian particle distributions, the running diffusion coefficient (t) is expected to converge in time to a constant value (t) → . Simulations are stopped after several orders of gyrations once the running diffusion coefficients converge, and the final value of the diffusion coefficient is taken as the steady-state one. We repeat these simulations several times with 2000 particles and the same parameters but varying random numbers for turbulence generation. As the low turbulence levels in the RSR require more statistics, we use 50 random seeds in the RSR for b∕B ≤ 2 and otherwise 20 random seeds. We perform this step for different reduced rigidities and fit a (8) power law ∝ to the RSR, where particles experience resonant scatterings: l min ∕( l c (b∕B)) < < l max ∕(2 l c ) . We determine the reduced rigidity scaling for a broad range of turbulence levels 0.05 ≲ b∕B ≲ 10 , while B = 1 G remains constant.

Simulation results
Over the last decades there have been great advances in direct numerical calculations of particle transport in mean magnetic fields subject to turbulent perturbations with a Kolmogorov-like spectrum accompanied by the adaptation of theoretical models to match these simulation results. In addition to studies of parallel diffusion coefficients, the perpendicular diffusion coefficients and the relationship between the two components have been the subject of research [12, 15-17, 28-32, 35, 38-43].
Good agreement between analytical models such as QLT (in some regimes), nonlinear guiding center theory and unified nonlinear theory with test particle simulations was found for two-dimensional, slab, or composite (twodimensional & slab) turbulence (see e.g. [44] for a review). Despite these advances, for isotropic three-dimensional turbulence, tension with these theories was found at small reduced rigidities [12,32].
In Fig. 2 we show simulation results of the parallel (solid lines) and perpendicular (dotted lines) diffusion coefficients calculated using 2000 particles in each simulation. For each of these ratios, 20 different energies are simulated. The statistical uncertainties of the diffusion coefficients are calculated by repeating each simulation 20-50 times with different realizations of the turbulent magnetic field. The averaged diffusion coefficients are shown as functions of the reduced rigidity . The statistical uncertainties are too small to be visible.
The classification of the diffusive transport into different transport regimes based on the particle energy and the turbulence level, and especially the definition of the resonant scattering regime (RSR), allows the comparison of simulated data with the QLT [12]. This is due to the fact that only in the RSR the conditions of QLT for resonant scattering over the whole pitch angle range are met in the simulations. As the lower RSR boundary = l min ∕( l c b∕B) utilizes the approximation of small turbulence levels, we apply this formula to determine the lower RSR boundary of the lowest used turbulence level b∕B = 0.067 . Since this lower boundary provides the smallest possible error from the underlying approximation for small turbulence levels, we also use the corresponding energy E = 7.5 PeV for all other turbulence levels as the lower boundary. While this approach does not exploit the full width of the RSR, as the RSR extends to smaller reduced rigidities as the turbulence level increases, the approximation of small angles still applies to all fits. The small errors in the fits are an indication that our chosen range in reduced rigidity is sufficiently large. The upper RSR boundary is = 5∕(2 ) for all turbulence levels.
Note that the simulated diffusion coefficients shown in Fig. 2 show qualitative similarities with comparable studies for isotropic three-dimensional turbulence (see e.g. [17]), but also with two-dimensional turbulence for the perpendicular component at large reduced rigidities [45].
In the following, first the parallel and then the perpendicular diffusion coefficients are investigated in the RSR and the QBR.

Parallel component
Especially for lower turbulence levels, the parallel diffusion coefficients are several orders of magnitude larger than the perpendicular components, as presented in Fig. 2. This is due to the different scaling with the turbulence level. While the parallel component scales with (b∕B) −2 , the perpendicular component decreases with decreasing turbulence level. Analytical theories describe the scaling of the turbulence dependence of the diffusion coefficients (see e.g. [26,33]). We perform fits in the RSR and QBR to determine the energy scaling of the diffusion coefficients for different ratios of b/B, as theory predicts a proportionality between log and log in both regimes. Note that due to the limitations of the fluctuation range of the synthetic turbulence, we can perform the fits of the predicted power law in the RSR only over a relatively small energy range, however, while strictly following our derived boundaries of this energy regime.
The calculated power-law indices ∥ of the fits ∥ ∝ ∥ are presented for different turbulence levels in Fig. 3. The simulation results for the QBR are shown as red crosses. The ∥ values are in agreement with theoretical predictions of ∥ = 2 [15,46].
For the RSR, we compare our results with simulation results from [12]. The underlying simulation setups of the shown simulation results from [12] and this publication differ only in the spacing of the grid points for the generation of turbulence. While in [12], a grid spacing of l min ∕10 was chosen, here we use the maximum possible grid spacing of l min ∕2 , which is based on the sampling theorem. The grid size has a direct influence on two numerical effects: Smaller spacing between grid points improves the resolution of small scales and thus reduces the numerical errors when interpolating the magnetic field between surrounding grid points. And, a larger grid reduces the need to periodically continue the original box when particles diffuse spatially across the boundaries. A comparison with a gridless method for turbulence generation has provisionally shown good agreement with the larger grid spacing [47]. The gridless method is based on [16,36,48] and evaluates at each point in space the sum of pre-generated wavemodes while eliminating the need of storing discrete field vectors on the grid.
While the results of both setups are slightly shifted against each other, they follow the same trend: The QLT limit of = 1∕3 is not yet reached at the 5% turbulence level, but is expected at even lower b/B based on the trend. At b∕B ≈ 2 , the upper limit for diffusion is reached, consistent with Bohm diffusion ∥ ∝ .
Here we would like to mention that for the fits at high turbulence levels an increasing number of simulated diffusion coefficients at low rigidities are considered, as shown in Fig. 2. Towards lower values and thus smaller scales on which the particles gyrate, numerical effects such as interpolation errors increase, which is due to the limited resolution of the discrete grid on which the turbulence is stored, which is why the values towards high turbulence levels may contain an additional systematic error, which is not indicated, but may lead to the deviation from the expected ≈ 1.

Perpendicular component
Fits based on ⟂ ∝ ⟂ for the perpendicular components are shown in Fig. 4. For comparison, simulated values from literature are included. The slopes are quite sensitive to the range of considered as pointed out in [12,17], which is why the fitted slopes vary based on different fit ranges in different publications. Nevertheless, the same trend can be seen across all publications: Increasing turbulence levels lead to a coincidence of the parallel and perpendicular components, as ∥ and ⟂ converge towards 1 2 . This is expected since the directional component of the background field loses influence with increasing turbulence levels and diffusion becomes isotropic. As with parallel diffusion, Bohm diffusion must also occur for perpendicular diffusion with regards to the background field for large turbulence levels (Fig. 5). At lower turbulence levels, however, the increase in field-line random walk (FLRW) is expected, since the diffusion coefficient of the field lines scales proportionally to (b∕B) 2 . Since FLRW is only independent of energy when particles follow the field lines, there is a transition in behavior at roughly the energy at which the field lines separate by more than one gyroradius over the distance the particles travel in one gyroorbit. The observable decrease of the energy dependence of particle diffusion with decreasing turbulence levels, which is associated with a smaller value of ⟂ , can thus be explained by the increasing importance of particles diffusing through FLRW.

Discussion in the astrophysical context
What do our results imply in the context of astrophysics? In general, we can show that the diffusion tensor strongly depends on the turbulence level -a fact that is typically not taken into account in astrophysical simulations. A prominent example where this b/B-dependence can play a role is the high-energy signatures for our own Galaxy, where the Fermi satellite has measured gammaray emission at GeV energies [50,51]. It was possible to derive the contribution from hadrons and thereby quantify both the cosmic-ray proton number density and the index of the cosmic-ray proton spectrum as a function of galactocentric radius. These results are presented in panel a) of Fig. 8. From these data follows what is known as the cosmic-ray gradient problem. In particular, there is a gradient in the proton index, ranging from relatively flat spectra in the innermost ∼ 5 pc of the Galaxy ( E −2.3 − E −2.5 ) to a steep index of close to E −3.0 at the outermost radii. Recent interpretations of this feature include phenomenological models of a change in the diffusion index based on geometric effects of anisotropic diffusion [52], a galactocentric radius dependent diffusion tensor [5,[53][54][55], and nonlinear cosmic ray transport with scattering and advection off self-generated turbulence in combination with galactocentric dependent cosmic ray source distributions [56].
With the turbulence-dependent energy dependence of diffusion presented here, we can expand the anisotropic diffusion models presented in the literature. To simplify the transport model even further we neglect momentum diffusion pp = 0 and adiabatic effects ∇ ⋅ = 0 and assume a stationary state (see Eq. 1) As the number of individual sources that contribute to S change with the position in the Galaxy, so does S. This does, however, not change the spectral index systematically, as in diffusive shock acceleration, the spectral index mainly depends on the strength of the shock. Only the intensity of the signal changes as the source density varies with galactocentric radius. The diffusion term can be approximated by using the effective escape distance d ∥ and d ⟂ in the parallel and perpendicular directions, respectively, as the spatial dependence of the term Here, the factors have been identified by escape times This is a very rough simplification, neglecting the full radial dependence of the diffusion tensor. We are doing this as, at this point, we simply want to discuss the effect caused by the perpendicular and parallel components and interpret the escape times as radially dependent for now. A quantitative model will have to take into account the full spatial structure of these terms. In our model, the central question is: which of these two time scales will dominate at what radius? The diffusion time scale is given by the shorter time scale, For the two individual diffusion time scales, these dependencies are as follows: -Energy and b / B : As discussed in the earlier sections, the energy dependence of the two coefficients is Here, the dependence is (b∕B) ±2 with a negative and positive turbulence level index for parallel and perpendicular diffusion, respectively. Thus, as the time scales behave as −1 i , we have the following dependences for the parallel and perpendicular time scales: -Field direction/galactocentric distance: The escape distance for diffusive transport depends on the field direction, which changes with galactocentric distance: for parallel transport, the escape direction is along the field lines. Figure 6 shows the magnetic field in the Galactic center region (central 2 kpc with a height of 300 pc). The field is a combination of the global field first presented by [57], here used in a modified version [58], plus a component in the galactic plane presented in [59]. This combination of fields is necessary, whereas global field models omit the in-plane component in the center. It can be seen here that even with the in-plane component, the field lines are essentially pointing perpendicular to the Galactic plane in the Galactic center region. This can also be seen in Fig. 7, where the mean angle of the field lines at a given galactocentric radius with respect to the plane direction is shown. Up to radii around 3 − 5 kpc, the angle is very close to 90 • , while at values larger than r gc ∼ 5 kpc, it becomes significantly smaller than 45 • . This is consistent with the presence of a Galactic wind in the inner 3 − 5 kpc [60,61], which can contribute to orient the field lines in the direction of the wind speed, which is perpendicular to the plane. In the outer region at r gc > 5 kpc, escape via parallel diffusion therefore preferentially occurs along the plane. The relevant escape distance for parallel transport is therefore a function of galactocentric radius. As a simplification, the parallel escape distance, defined as the length of the escape path in parallel along the mean magnetic field, can be approximated as the scale height d ∥ ≈ H ≈ 300 pc in the inner Galaxy, i.e. for r gc ≲ 5 kpc. In the outer Galaxy, r gc ≳ 5 kpc, it can be approximated as the inplane propagation distance. The latter is significantly longer, d ∥ > r max − r gc ≫ 300 pc for r gc > 5 kpc, with r max ∼ 20 kpc, especially when considering that particles will not diffuse out in straight lines along r , but rather follow the field lines. The perpendicular transport shows the opposite behavior: in the inner Galaxy, perpendicular to the field lines means in-plane Fig. 6 Three-dimensional view on the magnetic field lines in the Galactic center, using a combination of the global field [58] and an in-plane Galactic center component [59]  angle to plane [rad] only regular field,x ± σ Kleimann+(2019) with CMZ Field Fig. 7 Mean angle of the field lines with respect to the Galactic plane as a function of the galactocentric radius. The blue line shows the pure global field [58], the black line shows the combination of [58] with the Galactic Center field model of [59] propagation and thus d ⟂ > r max − r gc ≫ 300 pc. In the outer Galaxy, perpendicular transport approaches an orientation perpendicular to the plane, with d ⟂ ≈ H ≈ 300 pc. A larger escape distance is accompanied by a larger escape time, which leads to the fact that in the center, parallel transport results in the fastest escape of particles. In the outer Galaxy, perpendicular transport is accompanied by shorter escape distances until the escape distance along the plane in the outermost part of the Galaxy becomes comparable to the scale height. Thus, at the outskirts of the Galaxy, parallel transport will dominate again.
Combining these two effects, we get the diffusion time scale Here, the notation ⟨...⟩ refers to averaging over all field lines at a given galactocentric radius. It should be noted that d ∥ (r) and d ⟂ (r) represent functions at a fixed position in the Galaxy and averaging over the scale height z ∈ [−H, +H] needs to be performed for detailed results. We refrain from doing so at this point because we would like to keep the argument simple and these details are not expected to change the interpretation.
Note that inhomogeneities in B at the kpc scale are apparent in the Galaxy. This implies that drifts will occur with characteristic escape time scales on the order of drif t ∼ 10 13 ⋅ (E∕100 GeV) −1 yr based on drift velocities [62] along the z-direction within the setup of our toy model. While a complete quantitative picture of cosmicray diffusion must necessarily incorporate such drift physics, the conclusions at which the present work has arrived can still be expected to apply, given that these drifts cause bulk rather than diffusive motion and their time scales are many orders larger than those of diffusion and advection derived in the following.
We can now design the following toy model to investigate the dependence on the galactocentric radius: based on the parameters that are displayed in Fig. 8, there is a change in the physics of the turbulence at r 0 ∼ 5 kpc that should be reflected in the equations. Here, we use a simple ansatz based on the turbulence level shown in Fig. 8: We can also parametrize the total magnetic field in a similar way, Here, we roughly approximate the total magnetic field behavior with r gc based on the function shown in Fig. 8 and use ∼ 0.1 in the inner 5 kpc and ∼ 1 for r gc > 5 kpc.
For the diffusion indices, we use constant values, as the dependency is rather small for the individual indices, We approximate the escape distances with the scale height for r < 19 kpc, since that is the dominant escape direction. For larger radii, d ∥ depends on the position in the Galaxy.
This result implies a somewhat steeper energy spectrum in the inner Galaxy, n ∝ E − s − i ∼ E −2.9 , using s ∼ 2.2 − 2.4 and ∥ ∼ 0.7 as compared to the outer Galaxy, where perpendicular transport dominates ( ⟂ ∼ 0.4 ) and thus n ∝ E −2.7±0.1 . needed. Here, we want to emphasize that the measurements of diffusive emission in the Galactic center (see the two innermost data points in Fig. 8) are dominated by contributions from point sources and are further confounded by interpolation effects during the analysis process [50], which limits the constraints imposed by these measurements in the Galactic center. However, as the observed change in the spectrum is from a rather flat one in the inner Galaxy, n ∝ E −2.3 , to a steeper index in the outer region, n ∝ E −2.8±0.1 , other effects need to play a role in the inner Galaxy, while the steepening in the outer Galaxy could be due to a change to diffusive behavior in general. Even a change from perpendicular to parallel escape in the outermost part of the Galaxy could be relevant.
In the inner part of the Galaxy, a wind could play an important role as already discussed in [5,60,61]. We also approximate the convective term in terms of a loss time scale, via introducing the convection time scale It was shown in [60,61] that a wind with speeds around 500-700 km/s is present in the Galaxy. Using an approximate value of u ∼ 5 ⋅ 10 7 cm/s and a scale height for the wind of d conv ∼ 300 pc, we obtain a time scale of independent of the particle energy.
Here, the terminal wind velocity is used, which might be overestimating the actual wind velocity below z < 300 pc. A dedicated wind model based on a stratified disk instead of a flat plane would improve the results.
We can compare this time scale with the diffusive one in the inner Galaxy, Here, we use the scaling with energy from our simulation findings (see Eqs. (17), (18), and surrounding text). However, the simulated absolute diffusion coefficients of the particles from Fig. 2 would have to be extrapolated over several orders of magnitude to apply to our model according to the above formula. Note that the exact scaling depends on the cosmic ray energy and the galactocentric dependent magnetic field properties such as the field strenght and the correlation length of the turbulence. Due to our uncertainties for from the fits, the errors in the scaling over several orders of magnitude would become too large. Therefore, we normalize the time scale dif f (r 0 = 5 kpc, E = 10 GeV) ∼ 10 7 yr by using dif f (r Earth = 8.5 kpc, E = 100 GeV) ∼ 5 × 10 6 years as the value of the escape time scale as measured at Earth [1]. Using ∼ 0.1 for the inner Galaxy and ∥ ∼ 0.7 , results in 2 − ∥ = 1.93 . This way, we get The Galactic wind dominates as a loss process if conv < dif f , which is the case for energies This equation is fulfilled for cosmic-ray energies E ≲ 700 GeV at r 0 , and still for energies E ≲ 10 GeV in the inner central molecular zone (CMZ), i.e. the inner 200 pc. Even by taking the uncertainties of the actual wind speed and the normalization of the diffusion time scale into account, cosmic-ray transport will be influenced by Galactic wind in the CMZ. Therefore, the measured spectral index in the energy range measured with Fermi as shown in Fig. 8 can be dominated by the wind as well. Thus, a (22) conv ∼ 5 × 10 5 years, model in which convection dominates in the inner Galaxy, changing to diffusive escape along the perpendicular direction of the field around r gc ∼ 5 kpc, fits the data well, with a possible further change to parallel transport in the outermost Galaxy. One prediction of this model is that the transport in the inner Galaxy should change from convection to parallel diffusion for energies larger than ∼ 100 GeV. This implies a significant change in the energy behavior of the diffuse spectrum from E −2.3 to E −3.0 at TeV energies. With the emerging Cherenkov Telescope Array (CTA) [63] and its use of more than a hundred IACTs in the Northern and Southern Hemispheres, the background will be reduced and the field of view increased to ≈ 10 • , while improving sensitivity by an order of magnitude compared to current observatories consisting of only a few groundbased Imaging Air Cherenkov Telescopes (IACTs) [64]. These enhancements from diffuse measurements of the inner Galaxy with CTA, once available, will allow verification of our model. Complementary high-sensitivity observations of Large High Altitude Air Shower Observatory (LHASSO) [65] and of the Southern Wide-field Gamma-ray Observatory (SWGO) [66] at higher energies will provide further insights.

Summary and conclusion
In this paper, we perform simulations of cosmic-ray diffusion for a range of particle energies and turbulence ratios b/B. We analyze the simulation results by only using data that fully lie in the resonant scattering regime, i.e. where diffusion can develop fully. In doing so, we find that a power-law fit performs well for the parallel and perpendicular components of the diffusion coefficient, i.e. i ∝ E i , with i = ∥, ⟂ . We also find that these spectral indices are functions of the turbulence level, i = i (b∕B) . In both cases, low turbulence levels lead to the lowest indices, ∥ (b∕B = 0.07) = 0.57 and ⟂ (b∕B = 0.07) = 0.33 . Both indices asymptote to i = 1 for values b∕B > 1 . This limit is consistent with expectations for b∕B → ∞ , as Bohm's theory predicts an index of 1 for purely turbulent fields. We show that the limit of QLT, which predicts an index for parallel transport of ∥ = 1∕3 , is not reached at turbulence levels b∕B ≳ 0.07 . Since the parallel diffusion coefficients increase significantly as the turbulence level decreases, simulations with even lower turbulence levels require a considerable amount of time for the plateaus of the running diffusion coefficients to develop. Potential numerical effects could then gain influence, which is why much smaller turbulence levels were not considered further [12]. The gradient of our result indicates that an index of ∥ = 1∕3 will be reached at much lower turbulence levels. There, the influence of field-line random walk on the effective perpendicular diffusion process increases. Since the FLRW is energy-independent, ⟂ also becomes smaller for decreasing turbulence levels.
Finally, we apply these results from fundamental plasma physics to the Galaxy. Here, measurements indicate a steepening of the local cosmic-ray energy spectra along the galactocentric radius, with The Kolmogorov approximation in QLT of a diffusionbased change in the energy behavior corresponding to a steepening by E −1∕3 would certainly be compatible with the inner, flat spectra. However, we were able to show here that for the turbulence level in the inner Galaxy, diffusion would steepen the spectrum significantly more, e.g. by E −2∕3 . Our conclusion is therefore that transport in the inner Galaxy must be influenced by a Galactic wind that does not change the energy spectrum, so that the spectrum after transport corresponds to the one at the sources for r gc < 5 kpc, i.e. n ∝ S ∝ E −2.35±0.05 . The large uncertainties of gamma-ray emission measurements in the inner Galaxy currently provides flexibility in the interpretation of these data, as evidenced by the interpretation of [56] of a very soft spectrum and of [55] and the current work of a hard proton spectrum in the inner part of the Galaxy. To resolve this tension, better data of the observationally challenging inner Galaxy is needed.
The steepening of the spectral index toward the outer Galaxy would then be due to the change in the transport process, from convective to diffusive. Here, we show that in the next-outer region, diffusive escape is dominated by the perpendicular transport component, with the perpendicular diffusion coefficient steepening the spectrum by E −0.4 , leading to a spectrum after transport of n ∝ S −1 ⟂ ∝ E −2.7 . A dominating perpendicular transport for these galactocentric radii is also suggested by e.g. [52,55]. At the outermost radii of the Galaxy, at around 20 kpc, diffusive escape will again move toward parallel transport and, the spectrum would steepen further. We show in this paper that these findings are consistent with observations at GeV energies. In contrast to [55], we expect a change in the spectrum at around 100 GeV-1 TeV cosmic-ray energies in the inner Galaxy due to the replacement of convective transport by parallel diffusion as the dominant escape process, which should result in a change from n ∝ S ∝ E −2.35±0.05 to n ∝ E ≈−3 according to our simulation results. In the future, more detailed simulations, accounting also for adiabatic energy changes that are currently neglected, and improved observations of the diffuse gamma-ray component in the Galaxy with CTA and SWGO over a greater energy range can help to test this picture and will help to discriminate the three discussed models.
Funding This work is supported by the "ADI 2019" project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02 (PR). PR also acknowledges support by the German Academic Exchange Service and by the RUB Research School via the Project International funding. LM acknowledges financial support from the Austrian Science Fund (FWF) under grant agreement number I 4144-N27. EZ gratefully acknowledges support by NSF AST2007323.

Data availability
The methods, results, and data in this paper are fully available to interested researchers.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.