Anomalous diffusion of cosmic rays: A geometric approach

The transport of cosmic rays is a phenomenon that fascinated many famous scientists (Fermi,¹ Ginzburg,² and many others) who gave fundamental contributions to tackle its complexity.³ In most of the approaches available in the literature, the transport of cosmic rays in the interstellar space is modeled using the Brownian motion framework that, in the continuum limit, leads to the classical diffusion equation.^3,4 This modeling approach is justified under the hypothesis that cosmic rays originating in the interstellar space are confined due to the magnetic field action, preventing them from leaving the galaxy. In this representation, also known as the isotropic diffusion model,⁵ the cosmic rays are free to travel in a quasi-homogeneous medium (the interstellar space) and, after interacting with localized scattering fields (e.g., shock waves and fluctuating hydromagnetic fields), are accelerated and scattered in random directions.³ The evolution of the cosmic rays’ energy spectrum with the distance from their source following these scattering processes is described in this model via the diffusion coefficient, which is approximated by a power-law function.

Nowadays, there are a number of studies using a different approach than the majority of galactic cosmic ray transport models relying on the assumption of isotropic diffusion. In this respect, the authors of Ref. 6 developed a numerical scheme that focuses on adapting the galactic matter distribution via consideration on the radiation field⁷ and more complex distributions of cosmic ray sources.⁸ They also provided a first implementation of a source localization model based on the H.E.S.S. Galactic plane survey. Another approach to describe the anisotropy of cosmic ray transport was provided in Ref. 9, where the authors implemented an algorithm to solve the transport problem using stochastic differential equations. This approach allowed the use of a diffusion tensor anisotropic with respect to the magnetic background field. In addition, in Ref. 10, the anisotropy of cosmic ray diffusion due to the large-scale structured galactic magnetic field was a key aspect in developing an accurate transport model. In order to include this effect, together with continuous loss processes, the authors developed a quantitative model, based on stochastic differential equations, of galactic propagation for cosmic ray protons.

In Ref. 11, the authors remark that the diffusion of cosmic rays becomes anisotropic when, on scales comparable with the Larmor radius of the particles, the intensity of the turbulent field (δB) is much smaller than the mean large-scale magnetic field B₀ (i.e., if δB/B₀ ≪ 1). Under this condition, the particles mainly diffuse along the magnetic field lines¹² and, in the limit of vanishing perpendicular diffusion coefficient, the cosmic ray transport across the main field occurs as a consequence of the (random walk-like) drifting magnetic field lines.¹³ By means of non-dimensional analysis for an idealized case of particles escaping a supernova remnant and diffusing into a magnetic flux tube with long coherence length (i.e., preserved for long distance), the authors show that the discrepancy between the diffusion coefficient estimated under the (opposite) assumptions of isotropic and 1D diffusion is much larger than an order of magnitude. These observations clarify the importance of the assumptions on cosmic ray diffusion to correctly interpret gamma-ray observations. Consequently, the authors of Ref. 11 developed a model for cosmic ray propagation assuming that they are highly magnetized and diffuse uniquely along the (drifting) magnetic field lines. The behavior of the energetic particles under the concurrent effect of particle diffusion along the magnetic field lines and the random walk of the very field lines is also known as anomalous diffusion. A comprehensive review on the models of anomalous diffusion of cosmic rays and energetic particles can be found in Ref. 14. Some of these research areas include energy transport in magnetized plasmas,^15,16 transport of energetic particles in solar wind,^17,18 acceleration of cosmic rays at astrophysical shocks,¹⁹ and their confinement in the galaxy.^20,21

More recently, the study of the motion of cosmic rays in random magnetic fields leads to the appearance of fractional operators in the cosmic-ray phenomenology^3,5 and to their formal definition. The application of fractional derivatives in the study of complex transport phenomena gained a great momentum also because of their connection with the stable non-Gaussian distribution.⁵ The most important feature of these new solutions is their power-law asymptotic behavior, which agrees well with the known properties of the turbulent interstellar medium, the Fermi acceleration mechanism, and other processes affecting cosmic ray propagation.³ In this work, we propose a stochastic scattering model devoted to the study of the behavior of cosmic rays propagating through the interstellar space and interacting with localized magnetized clouds (i.e., self-generated magneto-hydrodynamic waves). The propagating cosmic rays travel undisturbed through the interstellar medium according to long collisionless trajectories (Lévy flights)^22,23 interrupted by chaotic interaction with the magnetized clouds, after that they are partly re-emitted in random directions and partly remain trapped. We intend to show that this intermittent mechanism leads to the occurrence of anomalous transport dynamics. Under such circumstances, to an observer placed on the observation plane in Fig. 1, the cosmic ray distribution would appear non-Gaussian with strong heavy-tails characterized by power-law behavior.

The capability to analyze the information content of the non-Gaussian (heavy-tailed) distribution could provide valuable insights into the nature of the complex and highly scattering environments crossed by the cosmic rays. This ability would in turn allow the development of effective remote sensing techniques that would greatly improve the acquisition of high resolution astronomical images.

In this section, we describe a simplified numerical model aimed to identify the physical mechanisms determining the occurrence of anomalous transport associated with cosmic ray propagation through the interstellar medium. The problem under investigation is schematically represented in Fig. 1. The model builds on that developed by the authors of Ref. 24 and is composed of two parts: (1) a 2D geometric Monte Carlo model that describes the transverse spreading of the cosmic rays as a function of the interstellar magnetic field density and (2) a 1D random walk model that, using parameters projected from the 2D model, makes it possible to identify the effect of non-local (long-range) interactions and of re-emission mechanisms occurring during the scattering process. We remark that the 1D model represents a cross section of the cosmic ray distribution disk on the projection plane A (see line L in Fig. 1).

This section reports the details of a stochastic scattering model intended to describe the transverse diffusion associated with cosmic ray propagation through the interstellar medium. For simplicity, the analysis is performed into a 2D domain (even though it can be easily generalized to the 3D case). The main objectives of the model are (1) to analyze the role of the interstellar magnetic cloud density in the transverse diffusion and (2) to generate the probability density functions to be used in the 1D random walk model for characterizing the long-range interactions. These probability density functions will be synthesized using a purely geometric approach based on the concept of “view factors.”²⁵ 

The physical problem considered in this work, namely, a control volume representative of the interstellar medium embedding (scattering) magnetic clouds, is schematically represented in Fig. 1. The scattering magnetic clouds are assumed to be randomly distributed and are modeled as tubes. This hypothesis reflects the dynamics of the magnetic lines of force that, because of self-entanglement, create magnetic tubes.³ The uniform background corresponds instead to the interstellar medium where the cosmic rays can freely propagate. We remark that the magnetic tubes are composed of self-generated and entangled magnetic lines, which cause the intercepted cosmic rays to partly trap while randomly re-ejecting the remainder. The stochastic model proposed in this work is mainly intended to capture the transverse diffusion occurring upon propagation of the cosmic rays through the medium. As already highlighted in Ref. 3, this is the principal mechanism leading to the formation of heavy-tails. For simplicity, the analysis will be performed in two dimensions, i.e., on the cross section A′ of the control volume.

In the model presented, the magnetic tubes are assumed to be randomly distributed and characterized by constant diameter d. The magnetic cloud density is calculated via the definition of surface fraction, i.e., ratio of tubes’ total surface over the total control surface (SF = N/4 · (d/H)², where N is the number of magnetic tubes and H is the radius of the control surface). In the proposed model, the initial interaction of the cosmic ray beam with the magnetic clouds is approximated by assuming the origin of the rays in the central tube. The model is intended to simulate the multiple scattering of the rays’ beam due to its interaction with magnetic clouds and their resultant transverse diffusion that leads to the occurrence of heavy tailed distributions.²⁴ 

As anticipated in Sec. II A, the developed model simulates the interaction of a cosmic ray beam with magnetic clouds according to a simplified stochastic ray tracing model. The number density of magnetic clouds is N, and they are randomly distributed over a circular control surface of radius H. The tubes’ diameter and the radius of the circular control surface are set to d = 1 and H = 100, respectively. The specific value of the ratio d/H in the following simulations was selected to ensure the permanence of the scattered rays within the control domain and to avoid boundary effects. The view factors associated with each magnetic tube were calculated analytically with the formula provided in Ref. 25. This model allows us to calculate analytically (i.e., using a purely geometric approach) the relative scattering between each couple of tubes in the domain as schematically shown in Fig. 2. The analytical formula used to calculate the view factors between couples of cylinders is derived in Ref. 25 and reported as follows:

This formula can be further simplified under the assumption of magnetic tubes with constant diameter d,

Several examples of 2D view factors distributions are shown in Fig. 3. These distributions illustrate the relative scattering interactions between each couple of magnetic tubes in the control domain. In particular, Fig. 3 shows the 2D view factors distributions for three different values of the surface fraction SF₁ = 1.25% (a), SF₂ = 5% (b), and SF₃ = 31.25% (c). Results are provided in a three-dimensional view. The vertical axis reports the magnitude of the relative view factor between each couple of tubes calculated with the analytical formula provided in Ref. 25. These view factors will be used to synthesize the probability density functions of the 1D model developed in Sec. III. Figures 3(b), 3(d), and 3(f) report the view factors of the central magnetic tube with respect to the remaining tubes for the cases of tube diameters 1, 2, and 5, respectively.

As discussed so far, the main aim of this work is the analysis of the anomalous diffusion of cosmic rays when propagating through randomly distributed magnetic clouds. For this purpose, this section is devoted to projecting the bi-dimensional model of magnetic cloud distribution to an equivalent 1D model. The underlying hypothesis to this projection operation is the invariance of transport characteristics for any cross section. The black dashed line L in Fig. 1 shows the equivalent 1D model. The underlying mechanism that controls rays’ scattering through magnetic clouds is unaffected by this dimensionality reduction.²⁴ 

In order to introduce the equivalent 1D model, let us consider the following random walk problem, where N particles are distributed over M boxes aligned in a single dimension (see Fig. 4). No limitation is assumed for the number of particles contained in each single box, so no limitations are imposed on the local particle density. Note that the particles in each box represent the rays forming the initial cosmic ray beam while the boxes represent the scattering magnetic clouds. As the particles encounter material boundaries, they are scattered in different directions, hence in different boxes. We are interested in the transverse scattering produced when the rays interact with the scattering clouds. In fact, it is this transverse scattering that induces the beam spreading and the transition to a non-Gaussian intensity distribution.

At each scattering event, the rays can jump into any of the remaining boxes according to a given probability distribution P(k − i), which is a function of the relative position between boxes. P(k − i) indicates the probability that a particle, initially in the i-th box, will jump into the k-th box after a single scattering event. The number of particles f in the box i between two consecutive scattering events n and n + 1 can be expressed as

Note that P must also satisfy the normalization condition

The four sums on the right-hand side of Eq. (3) represent, respectively: (1) the particles entering the box i from the left, (2) the particles entering the box i from the right, (3) the particles escaping the box i to the left, and (4) the particles escaping the box i to the right. The assumption that the number of particles is much larger than the number of boxes has two important consequences: (1) it reduces considerably rounding errors and (2) it allows a discrete problem to be treated with a continuum approach. The non-dimensional constant A in Eq. (3) can be interpreted as the total fraction of particles that is scattered and thus is subject to the condition A ≤ 1. For example, at any given time if A = 1, all the particles in any box are removed and distributed among the other boxes, if A = 0.5 only half of the available particles is removed. The coefficient A could also be dependent on either the scattering event or the specific box location (i.e., $A=Ain$⁠) in order to account for inhomogeneity of the medium.²⁴ At this stage, we consider A as a constant. Note also that this non-dimensional parameter could potentially be a function of the number of particles. However, such dependence would make the governing equation non-linear therefore requiring a different approach.

In the following analysis, we will consider three different values of A, i.e., 0.1, 0.5, and 0.9 corresponding to magnetic clouds characterized by an increasing probability to re-emit the intercepted cosmic rays. Assuming that P is symmetric, i.e., if P(j) = P(−j), the balance equation can be recast as

which implies that no net exchange occurs between two boxes containing the same number of particles $(fkn=fin)$⁠. A further manipulation of Eq. (5) yields

where M_L and M_R are, respectively, the number of boxes on the left and on the right of the i-th box. If the number of boxes is very large (M → ∞), assuming that the distribution P(j) → 0 for j → ±∞, and invoking the normalization condition Eq. (4), we obtain

which simply states that the number of particles in the box i after the scattering event is equal to the number of particles already in i and that are not scattered away, plus the particles scattered into i from other boxes.

The model introduced in Sec. III is devoted to capturing the characteristics of the transverse scattering of cosmic rays interacting with magnetic clouds. Each box of the 1D model (that we recall is a projection of the 2D configuration) represents a magnetic cloud. In this work, we considered random configurations of scattering cylinders and we introduced concentric levels to perform the integration of the view factor distribution. This dimensional reduction corresponds to the assumed isotropic cross section. The probability density function (pdf) of a given box is therefore given by the average of the view factors of a given level of the 2D original configuration.²⁴ 

In addition, in the equivalent 1D model (as in the original 2D configuration), all the rays are contained initially in the center box at the initial time instant. The interaction between boxes at successive time instants is determined by the view factors (i.e., the pdfs).

We performed a parametric study where we analyzed the evolution of the ray distribution (injected in the central box) as a function of the surface fraction. We observe that with increasing surface fraction (that is, increasing the diameter of the scattering cylinders), the rays’ distributions become α-stable (non-Gaussian) with a strong superdiffusive character, as illustrated in Figs. 5–7. We remark that in these figures, the red circles correspond to the calculated data while the dotted blue lines are the fit to an α-stable distribution. Plots are presented on different scales in order to highlight the behavior of the distributions in both the central and tail regions. These results clearly show the existence of ranges of density and re-emission strengths of the magnetic clouds that result in heavy tailed distribution of the propagating cosmic rays. This denotes the occurrence of transport dominated by long-range interaction mechanisms. We recall that in the 1D random walk model, the occurrence of long-range interaction mechanisms corresponds to rays able to jump several boxes in one step, and at a macroscopic level this results in anomalous diffusion. The values of the parameter α obtained from the α-stable fit of the numerically generated distributions for three different values of the surface fraction are summarized in Table I.

We remark that the values and the global trend of the stability index α obtained in the previous numerical simulations closely resemble the trend of the power-law scaling exponent reported in Ref. 3.

We recall that the parameter α is the characteristic exponent of an α-stable distribution and it defines its degree of impulsiveness. A decreasing value of α corresponds to a distribution with increasing value of impulsiveness.²⁶ The main feature of an impulsive distribution is the presence of marked heavy-tails when compared to a Gaussian distribution (details in the Appendix of Ref. 27). The use of α-stable distributions (and the related statistics)^26,28,29 has already proved a powerful tool in analyzing astrophysical data with the aim of characterizing extragalactic sources. The approach proved valid in broad frequency ranges (including microwaves,^30–32 x rays, and radio-astronomy³³) and showed its robustness when dealing with Gaussian-noise contaminated data.²⁶ In this perspective, the statistical mechanics approach proposed in this work provides a solid framework to connect the parameters of the α-stable distribution to the spectrum of the galactic structure properties including density and composition of the magnetic clouds, and the intensity and location of cosmic rays sources. The ability to decode this information is key in many astrophysical scenarios.

Besides, the α-stable heavy-tailed distributions observed in this work, by means of the connections they have with fractional order governing equations,^24,27,34 open the possibility to model intermittency and therefore to capture the dynamical characteristics of highly complex structures like the magnetic clouds present in the interstellar space.

It is well known that one of the first applications of the Monte Carlo method is due to Enrico Fermi while studying neutron diffusion. His calculations based on Monte Carlo integration would eventually lead to the Quantum Monte Carlo method (i.e., diffusion Monte Carlo methods). In this perspective, the results of the geometric Monte Carlo technique used in this work point to a novel approach of performing integration in higher dimensional spaces. Several works^3,27 have already pointed out the correspondence existing between the parameters of the α-stable fit and Green’s function of the governing fractional Cauchy problem. In this section, we perform the fit of the distribution reported in Fig. 5(e) (i.e., SF = 1.25% and scattering probability A = 0.9) using Green’s function of the governing fractional diffusion equation, as discussed in detail in Sec. 7 of Ref. 27. In the approach proposed in Ref. 27, the authors exploited the connection between the parameters of the α-stable fit and the governing fractional Green’s function to obtain a closed-form solution to the decay envelope of the acoustic intensity. The fractional Green’s function (reported in Fig. 8) can then be easily integrated (the corresponding integral area is highlighted in Fig. 8 with vertical blue lines bounded below by the x axis and above by Green’s function). We remark that this integration technique contains all the ingredients that are typically necessary for Monte Carlo integration, which are: (1) a source, (2) scattering objects, and (3) properties of the scattering objects. On the other hand, this technique can exploit the connection between the parameters of the α-stable fit and the fractional Green’s function (that is a closed-form analytical solution) to provide a new efficient tool for performing integration in higher dimensional spaces.

This work proposed a simplified statistical mechanics approach to analyze the anomalous diffusion of cosmic rays through randomly distributed magnetic clouds. The results of the geometric-based Monte Carlo simulations showed the emergence of anomalous diffusion for such transport processes. The characteristics of the anomalous transport were parametrically studied for wide ranges of magnetic cloud density and re-emission strength. The analysis revealed the ranges of parameters where cosmic ray distribution showed heavy-tailed nature, i.e., with α-stable characteristics. This marked deviation of the tail region from the Gaussian distribution highlights heavy information content that, once properly decoded, could provide valuable insights into the characteristics of the magnetic cloud structure. These results could be key in the challenge of achieving high resolution astrophysical images that require remote sensing technologies able to interpret the transport characteristics in the complex interstellar environments.

The authors thank Dr. Damiano Baccarella for comments that improved the manuscript.

The data that support the findings of this study are available from the corresponding author upon reasonable request.