# Cosmic Ray Small Scale Anisotropies and Local Turbulent Magnetic Fields

###### Abstract

Cosmic ray anisotropy has been observed in a wide energy range and at different angular scales by a variety of experiments over the past decade. However, no comprehensive or satisfactory explanation has been put forth to date. The arrival distribution of cosmic rays at Earth is the convolution of the distribution of their sources and of the effects of geometry and properties of the magnetic field through which particles propagate. It is generally believed that the anisotropy topology at the largest angular scale is adiabatically shaped by diffusion in the structured interstellar magnetic field. On the contrary, the medium- and small-scale angular structure could be an effect of nondiffusive propagation of cosmic rays in perturbed magnetic fields. In particular, a possible explanation of the observed small-scale anisotropy observed at TeV energy scale, may come from the effect of particle scattering in turbulent magnetized plasmas. We perform numerical integration of test particle trajectories in low- compressible magnetohydrodynamic turbulence to study how the cosmic rays arrival direction distribution is perturbed when they stream along the local turbulent magnetic field. We utilize Liouville’s theorem for obtaining the anisotropy at Earth and provide the theoretical framework for the application of the theorem in the specific case of cosmic ray arrival distribution. In this work, we discuss the effects on the anisotropy arising from propagation in this inhomogeneous and turbulent interstellar magnetic field.

###### Subject headings:

cosmic rays - magnetic fields - magnetohydrodynamics (MHD) - scattering - turbulence## 1. Introduction

Cosmic rays are found to possess a small but measurable anisotropy in their arrival direction distribution at Earth. The origin of the observed anisotropy is not yet understood. However, it is reasonable to assume that it is a combination of effects correlated to the distribution of the galactic sources of cosmic rays, the geometry and turbulence properties of the galactic magnetic field, and the propagation in interstellar magnetized plasmas. These are likely to be responsible for the complex shape of the energy spectrum as well (Gaisser et al., 2013). Since we don’t know the locations of cosmic ray sources and we lack details of the interstellar magnetic field, understanding these observations of anisotropy is not an easy task.

Above the energy range where cosmic rays are directly affected by inner heliospheric processes (see, e.g., Florinski et al. (2013); Manuel et al. (2014); Florinski et al. (2015)), a statistically significant anisotropy has been observed by a variety of experiments, sensitive to different energy ranges (from tens of GeV to a few PeV), located on or below the Earth’s surface in the Northern Hemisphere (Nagashima et al., 1988; Hall et al., 1999; Amenomori et al., 2005, 2006; Guillian et al., 2007; Abdo et al., 2009; Aglietta et al., 2009; Zhang et al., 2009; Munakata et al., 2010; Amenomori et al., 2011; de Jong et al., 2011; Shuwang et al., 2011; Bartoli et al., 2015) and in the Southern Hemisphere (Abbasi et al., 2010, 2011, 2012; Aartsen et al., 2013).

The global anisotropy changes with energy in a non-trivial fashion. From about 100 GeV to tens of TeV, it has been observed to have an approximately consistent structure at the largest scale, although its amplitude appears to increase with energy. Above a few tens of TeV, the observed progressive change in the anisotropy topology may indicate a transition between two processes shaping the particles’ arrival distribution at Earth. The observation could be qualitatively explained on the basis of diffusive propagation of cosmic rays in the Milky Way from stochastically distributed sources, responsible for generating a gradient in cosmic ray density.

Numerical studies of particle propagation in a scenario of homogeneous and isotropic diffusion in the Galaxy predicts that the cosmic ray density gradient, and the consequent induced anisotropy, has a dipole shape with direction toward the source of the particle and with amplitude that directly depends on the diffusion coefficient. In particular, for a given realization of galactic source spatial distribution, the dipole anisotropy would point toward the source with the largest contribution (Erlykin & Wolfendale, 2006; Blasi & Amato, 2012; Ptuskin, 2012; Pohl & Eichler, 2013; Sveshnikova et al., 2013; Savchenko et al., 2015), which may change with energy, in agreement with observations. On the other hand, the systematic overestimation of the anisotropy amplitude may be partially compensated by the fact that diffusion is expected to be anisotropic (see, e.g., Effenberger et al. (2012)), thus modifying the expected cosmic ray density gradient shape as a function of the source direction with respect to the regular galactic magnetic field (Kumar & Eichler, 2014). The misalignment between the cosmic ray density gradient and the regular galactic magnetic field would prevent pointing to any specific source, and it would suppress the anisotropy amplitude to a value closer to what has been observed (Mertsch & Funk, 2015).

Another scenario is that it is the transition from heliospheric- to interstellar-dominated contributions, starting at a 10 TeV energy scale and culminating around 100-200 TeV, at the origin of the shift in anisotropy. In Desiati & Lazarian (2013), this scenario was proposed noting that 10 TeV protons have an average gyroradius, in a G scale magnetic field, on the order of the thickness of the heliosphere (see, e.g., Pogorelov et al. (2006)). In addition, the dynamical instabilities of the heliospheric magnetized plasma at smaller scales may generate strong scattering that redistributes the arrival direction of TeV cosmic rays. This heliospheric scenario is studied and presented in our companion paper López-Barquero et al. (2016). Such strong scattering may be able to produce large localized particle gradients, experimentally interpreted as medium- or small-scale anisotropy. In Schwadron et al. (2014), a scenario of weak influence of the heliospheric magnetic field on TeV protons was explored, thus interpreting the observations as directly related to the ordering of the local interstellar magnetic field (see also Zhang et al. (2014)).

The anisotropy appears to possess a complex angular structure with evidence of a harder cosmic ray spectrum within the localized excess region in the apparent direction of the heliospheric tail (Amenomori et al., 2007; Abdo et al., 2008; Bartoli et al., 2013; Abeysekara et al., 2014). The decomposition of TeV cosmic ray anisotropy in the individual spherical harmonic contributions shows that the arrival direction distribution is dominated by large-scale structures (such as, e.g., dipole and quadrupole) with a relative intensity on the order of 10, but that medium- and small-scale angular structures are significantly contributing with relative intensities below 10.

The experimental determination of small angular scale anisotropy is generally performed by filtering out the large-scale modulations from the observed arrival direction distribution, thus retaining all structures with large angular gradients. Some of the small-scale anisotropy features seem to be correlated to regions in the sky where the global anisotropy has large variations (see Desiati & Lazarian (2013)) or may be an effect of re-acceleration by magnetic reconnection processes in the tail of the heliosphere (Lazarian & Desiati, 2010; Desiati & Lazarian, 2012). However, globally, the observed small scale anisotropy may appear to be rather randomly structured and, therefore, possibly a natural consequence of cosmic ray propagation in the local turbulent magnetic field in the presence of a global anisotropy (Giacinti & Sigl, 2012). The global anisotropy, at all angular scales, arises from the same physical processes, thus it is impossible to disentangle its origin. However, as a first approximation, it is generally assumed that the anisotropy at the largest scale is dominated by global physical processes (such as a density gradient from sources of cosmic rays or from convective effects originated by large-scale cumulative stellar winds, and from propagation through the regular galactic magnetic field), while the small angular scales are dominated by local processes.

A complex angular power spectrum is asymptotically generated by progressive decomposition of the energy of an initial anisotropy distribution (for instance a dipole) into higher multipoles, by the effect of scattering off magnetic turbulence (Ahlers, 2014; Ahlers & Mertsch, 2015). The conservation of phase space density, as stated by Liouville’s theorem, predicts, for an idealized situation of a homogenous large-scale anisotropy, the total sum of the multipoles’ angular terms is conserved. This makes it possible to generate small-scale structures, as shown in Ahlers (2014).

Turbulence in astrophysical plasmas have significant effects on particle propagation, in that the stochastic nature of magnetic field lines is transmitted to the particles’ trajectories. If particles are tied to magnetic field lines, the maximum perpendicular diffusion rate is set by the rate of perpendicular field line wandering, scaled by particle velocity (field line random walk). This has been extensively discussed by, e.g. Jokipii (1966); Jokipii & Parker (1969); Rechester & Rosenbluth (1978); Giacalone & Jokipii (1999); Minnie et al. (2009); Shalchi (2009) in the context of diffusion at distance scale larger than the turbulence injection scale. On scales smaller than the injection scale, particles are characterized by super-diffusion in the perpendicular direction of the mean magnetic field. This behavior is described by the so-called Richardson diffusion, where particle separation grows as (time) (Lazarian & Yan, 2014). In Eyink et al. (2011) it was shown that this is directly connected to the separation of turbulent magnetic field lines, which grows as (distance) (Lazarian & Vishniac, 1999). In this case, the stochastic nature of the astrophysical magnetic fields (such as, e.g. the interstellar magnetic field, or, at larger scales, the intercluster magnetic field), inevitably produces chaotic particle trajectories, meaning that the geometry of particle trajectories is highly sensitive to the actual initial conditions. Diffusion by magnetic field line wandering is found to be a dominant contribution and stronger than the extreme case of Bohm diffusion, where scattering frequency is one per particle gyration. The important aspect that determines the properties of a large ensemble of particles is its statistical nature, which in this case is influenced by the properties of the magnetic field and specifically by the induced scattering rate. While individual trajectories may have chaotic properties and are, therefore, practically/realistically unpredictable, the large ensemble of particles can still be deterministically described. For these reasons, the recorded spatial distribution of the ensemble will be statistically determined and a direct consequence of the properties of the turbulent magnetic field. In this work, a study of particle propagation in compressible magnetohydrodynamics (MHD) turbulence is performed in the context of its effects on arrival direction distribution.

In addition, depending on the degree to which magnetic field lines diverge on scales less than the particle gyroradius, pitch angle scattering on small-scale magnetic perturbations affects particle distribution at the large spatial scale, thus increasing the diffusion coefficient depending on the large geometrical scale of the magnetic field turbulence (see, e.g., Desiati & Zweibel (2014) and references therein). The anisotropic nature of interstellar turbulence and the properties of turbulence itself can significantly complicate the description of cosmic ray transport (see, e.g., Yan & Lazarian (2008, 2011); Beresnyak et al. (2011)).

This paper is organized as follows. In section 2, we describe the turbulent magnetic field used and the test particle trajectory integration. In section 3, we discuss the validity of applying Liouville’s theorem in the context of this work by statistically assessing the level of conservation of magnetic moment. The results of the study are presented in section 4 and discussed in section 5. Concluding remarks follow in section 6.

## 2. Cosmic ray propagation in turbulent magnetic fields

Cosmic rays, which are accelerated in their sources and “injected” into the interstellar medium (ISM), are free to propagate through the interstellar magnetic field. Globally, the galactic magnetic field is characterized by a large-scale regular component and a small-scale random component (see, e.g., Jansson & Farrar (2012a, b)). The regular component can be described as a superposition of spiral and toroidal structures, possibly with a contribution perpendicular to the galactic plane, and the random component represents the stochastic perturbations of the regular field caused by the dynamics of the galaxy and its density distribution. An important property of the interstellar magnetic field is turbulence. A variety of observations show that the coherence scale of turbulent magnetic fields is on the order of 10 pc in spiral arm regions (with more frequent stellar formation activity) and on the order of 100 pc in the interarm regions (Haverkorn et al., 2008). The injection scale of turbulence is determined by the scale at which the magnetic perturbation is generated (for instance, by stellar collapse or binary mergers).

Astrophysical plasmas are typically highly ionized and have high Reynolds numbers, thus the dynamics of the flow is dominated by nonlinear convective processes at the largest scale. In such conditions, turbulence develops and magnetized Alfvénic eddies dynamically cascade to smaller scales and progressively elongate along the magnetic field lines, as initially proposed by Sridhar & Goldreich (1994); Goldreich & Sridhar (1995) (see also Lazarian (2007) for a review). This model of incompressible MHD turbulence predicts a Kolmogorov-type energy power spectrum in terms of the wave-vector component perpendicular to the local direction of the magnetic field, while the parallel component of the wave vector is . In those pioneering papers, the theory assumes the injection of energy at scale and the injection velocity equal to the Alfvén velocity in the fluid , i.e., the Alfvén Mach number (i.e., trans-Alfvénic turbulence), where is the plasma velocity at injection scale . The model was later generalized for both sub-Alfvénic (i.e., ) and super-Alfvénic (i.e., ) cases (Lazarian & Vishniac, 1999; Lazarian, 2006) (see also Lazarian et al. (2012); Burkhart et al. (2014)). Typically, the ISM is characterized by 1 (Gaensler et al., 2011).This means that magnetic field lines do not typically fluctuate too far from the mean direction.

While the model by Goldreich & Sridhar (1995) describes incompressible MHD turbulence, modeling compressible turbulence turned out to be more complex. In isothermal plasmas, there are three types of MHD waves: Alfvén, slow, and fast waves. Alfvén modes are incompressible, while slow and fast modes are compressible. The compressible modes are conjectured to resemble incompressible behavior at high- (i.e., high gas to magnetic pressure ratio) (Lithwick & Goldreich, 2001); however, ISM plasmas are typically characterized by low-. Cho & Lazarian (2002) investigated the scaling properties of low- compressible sub-Alfvénic MHD turbulence and confirmed that slow compressible modes behave as the Alfvén incompressible modes but also that the fast modes are isotropic, since their velocity does not depend on magnetic field direction.

The spatial distribution of particles propagating in a turbulent field is affected by the magnetic perturbations within the particle autocorrelation length scale (or mean free path). To study the correlation between turbulence and cosmic ray distribution, trajectories of test particles were integrated in an MHD turbulent magnetic field. The methodology used is described in the next two sections.

### 2.1. Turbulent magnetic field

Test particle trajectories are integrated in a compressible sub-Alfvénic isothermal MHD turbulence in low- based on numerical calculations developed by Cho & Lazarian (2002). In the model, turbulence is driven solenoidally in Fourier space, setting the velocity and density fields initially to unity. The calculation is performed in a cube with side grid-points, with inertial range from grid-points (i.e., ) down to a “damping” scale of grid-points.

The average magnetic field is directed along the x-axis and turbulence is characterized by a gas-to-magnetic pressure value of and by Alfvénic Mach number . Such a Mach number corresponds to the fluctuations being on the order of the mean magnetic field at the injection scale, which is in agreement with what we would have expected from the local ISM. The external mean magnetic field is the only controlled parameter in this MHD model.

### 2.2. Cosmic ray propagation

The analysis is performed by integrating proton trajectories in a static magnetic field, using the set of 6-dimensional ordinary differential equations

(1) | |||||

(2) |

describing the Lorentz force and the particle velocity , where is the particle position vector and the momemtum. For , we use one steady state realization of the magnetic field described in section 2.1. The equations are integrated using the Bulirsch-Stoer integration method with adaptive time step. At each integration step, the magnetic field is interpolated using a 3D cubic spline, and integration is stopped when particles cross the border of the MHD simulation box. The accuracy of the particle trajectory integration for this study was assessed and is discussed in Appendix A. The choice of one specific realization of the magnetic field is justified by the fact that particle velocity is much larger than the plasma Alfvén velocity, thus induced electric fields can be neglected.

The magnetic field can be interpreted as a snapshot of the local interstellar magnetic field. Since the MHD simulation can be scaled to an arbitrary injection scale, in this study, two physical scenarios are investigated: an injection scale of = 10 pc, which is approximately the magnitude typical of the Milky Way spiral arms, and = 100 pc, which is characteristic of the interarm regions (see section 2). Although the ISM surrounding the solar system seems to have interarm properties (Frisch et al., 2012, 2014), both cases are taken into consideration. Assuming a particle gyroradius of 5 grid-points (corresponding to the smallest spatial scale in the MHD simulation), the injection scale sets the particle energy scale. In order to avoid grossly underestimating the effects of magnetic perturbations smaller than the particle gyroradius, the latter is set to 20 grid-points as well.

Set | Particles | |||||
---|---|---|---|---|---|---|

1 | 750 TeV | 0.24 pc | 10 pc | 5.0 pc | 3 G | |

2 | 7.5 PeV | 2.4 pc | 100 pc | 50 pc | 3 G | |

3 | 30 PeV | 10 pc | 100 pc | 60 pc | 3 G |

Table 1 shows the three trajectory integration data sets used in this study. They cover an energy range from 750 TeV to 30 PeV, partially overlapping the recent cosmic ray observations reported by the IceCube Observatory (Abbasi et al., 2012; Aartsen et al., 2013, 2015).

In order to study the effect of interstellar magnetic turbulence on the arrival direction distribution of cosmic rays, a large number of particles would need to be injected with randomly uniform directions on a spheric surface centered at the Earth and with a radius larger than the mean free path. Unfortunately, such method would be highly inefficient because a large fraction of the injected particles would never reach Earth. The alternative approach is typically to use the so-called “back-propagation method,” where the trajectories of particles with opposite charge are integrated from Earth, with initial directions uniformly distributed, backward into outer space. Since energy losses are negligible for proton particles, their energy is conserved, and therefore their trajectories can be time-reversed, provided there are no collisional scattering processes and no resonant scattering. Under such general conditions, Liouville’s theorem states that particle density in phase space is conserved along particle trajectories. In the presence of turbulence, magnetic fields can vary faster in space than the particle gyroradius, thus breaking adiabaticity and effectively inducing collisional processes. In section 3, the applicability of Liouville’s theorem in the present case is discussed.

For set 1 of Table 1, proton trajectories are integrated, while for sets 2 and 3, the number of particles is . The numerical trajectory integration starts from a point at the center of a 33 MHD simulation box, with uniform direction distribution, and stops when particles reach the edge of the integration volume. The integration box corresponds to a distance scale of 75 pc (for set 1 in table 1) and of 750 pc (for sets 2 and 3 in Table 1). The back-propagated trajectories calculated in this way provide the information on the effects of the magnetic field on an initially uniform particle distribution emanating from one point. In other words, they provide a map of regions that are magnetically connected to the origin at a given distance. The calculation implicitly takes into account the dynamical processes of a particle’s motion in a magnetic field, including drifts and pitch angle scattering. Under the hypothesis that such trajectories are time-reversible (see section 3), they can be interpreted as directly propagating from outer space back to the original point (assumed to coincide with Earth). This means that the particle distribution far from Earth, resulting from the back-propagation numerical calculations, corresponds to a perfectly uniform arrival distribution at Earth.

The numerically calculated trajectories can be used, therefore, to determine the arrival distribution at Earth as a consequence of a global anisotropy at a large distance. Such a global anisotropy, which is the effect of a particle density gradient induced by sources of cosmic rays or by convective processes, is treated as a weight on the forward-inverted trajectories, and it effectively produces a nonuniform arrival direction distribution at Earth, which is described in section 4.2.

## 3. The validity of Liouville’s theorem

Generating a large number of particle trajectories that pass through a point in space is implicitly highly inefficient, since most particles will bypass the target point. One possibility is to increase the size of the target to record those particles that pass “nearby,” or to appeal to Liouville’s theorem.

Liouville’s theorem states that the density of particles in the neighborhood of a given system in phase space is constant if restrictions are imposed on the system (Goldstein et al., 2002), as shown below. To determine its validity for our specific case of cosmic ray arrival distribution, we shall start with the continuity equation in phase space (Bradt, 2008a, b), under the assumption that the number of particles stays fixed:

(3) |

where is the del operator in momentum space, the distribution function, and the applied external force.

Upon expansion of this expression, we arrive at:

(4) |

The first term on the right-hand side of Eq. 4 vanishes, since the divergence of the velocity in phase space is zero. The third term in the right-hand side of Eq. 4, with , gives us restrictions on the forces that can be applied to the particles. For this term to go to zero, we need the forces to be conservative and differentiable. In particular, they must be p-divergence free. Evidently, if collisions are present, they will not fulfill these requirements. Now, we can analyze the case of magnetic forces and can express this term explicitly:

(5) |

Assuming that the magnetic field is independent of the velocity of the particle, the p-curl of goes to zero. In this case, the assumption is valid since the particles are moving so fast that the possibility of changing the magnetic field is negligible. For the p-curl of , it just vanishes, since the velocity is the gradient of a scalar function, the relativistic energy. Therefore, if we have only pure magnetic forces, the cancels.

We arrive at the equation:

(6) |

But this is just the expression for the total derivative of the distribution. Therefore, we can reexpress it as:

(7) |

which is the precisely the expression for Liouville’s theorem (Goldstein et al., 2002; Bradt, 2008b).

This derivation is for a pure magnetic force, but in fact when calculating the particles’ trajectories in a turbulent magnetic field, a variety of factors come into play. The most significant effect is when particles encounter a region where the magnetic field varies abruptly, i.e., the scale of variation of the magnetic field is shorter than the gyroradius of the particle. In this scenario, the trajectory does not have time to adjust smoothly to this change, and the interaction can be effectively considered a collision. For such cases, the right-hand side of Eq. 6 can be modified by the addition of a term, , which takes into account collisions of various origins that are differentiated by their exact functional form, given the fact that they will cause a nonzero time rate of change in the distribution function (Baumjohann et al., 1996). Therefore, under these conditions, Liouville’s theorem can’t be applied. To ensure that the abruptness in the trajectory is limited, we can calculate how adiabatic this change is. Having established this, it will ensure that a smooth variation will not greatly modify the density in phase space.

In the presence of collisions, the magnetic moment of the gyrating particles changes. Therefore, to check for the adiabaticity of the trajectories, we can calculate the magnetic moment for each particle at each time step and find out if it statistically behaves as an adiabatic invariant.

The relativistic magnetic moment (also called first adiabatic invariant) is given by:

(8) |

where is the momentum perpendicular to the magnetic field and the particle mass. This quantity, relating magnetic field and perpendicular momentum of the particle, is conserved if the field gradients are small within distances comparable to the particle gyroradius. Using the integrated particle trajectories from the data sets of Table 1, the magnetic moment Eq. 8 was calculated at each integration time step and histogrammed in time step slices. In each time slice, the mean value of the magnetic moment of the particle ensemble from each data set and the corresponding standard deviation were calculated. Figure 1 shows the evolution of over the integration time of the 30 PeV set (on the left) of Table 1 and the ratio / for the three sets (on the right) of Table 1. A perfect distribution with / = 0, indicates that the particles most likely stay in the same magnetic field line, which is unrealistic for a particle moving along a turbulent magnetic field. Nonetheless, a distribution peaked at a value much larger than one will tell us that the particles suffered strong variations in their trajectories and collision-like interactions happened. The distributions calculated for our three different energy data sets peak at around 0.5 and do not appear to have trends or large variations during integration time (compared to ), as shown on the left of Figure 1. This indicates that even though the particles have interacted with the turbulent field, and it has changed their trajectories, the changes are relatively smooth and with limited statistical impact on the overall particle ensemble.

Note that the width of the distributions in / means that at some level the magnetic moment of the particles fluctuates during propagation. Different particles follow different magnetic field lines and experience different interactions. There might be also a contribution from the limited accuracy of the numerical integration. However, as discussed and shown in appendix A, in the present application, the effects due to accuracy limitations are not sufficiently large to significantly violate adiabaticity. The highest level of possible inaccuracy, where particles are stochastically re-distributed at each gyration, provides a numerical diffusion at the level of Bohm diffusion, where particles are scattered every gyration. It has been proven that in astrophysical turbulence Bohm diffusion is always much smaller than diffusion induced by magnetic field line wandering (Lazarian & Yan, 2014). So the fact that the accuracy of the used trajectory integration method is significantly below the Bohm diffusion level (see appendix A) poses no problem on the statistical results obtained in this study.

Based on the fact that the ensemble-average / 1, especially that only magnetic forces are explicitly taken into account and that the changes in the particles’ trajectories are relatively smooth, it is assumed that Liouville’s theorem is applicable (see further discussions in Section 5). So in this study, back-propagation of the particle trajectories is justified in this statistical sense. The situation changes if scattering rate is higher, such as in the heliospheric magnetic field case studied in López-Barquero et al. (2016), where resonant scattering processes produce larger deviations from adiabaticity.

## 4. Results

This section shows the results obtained with the numerical calculation data sets described in Section 2.2.

### 4.1. Mean Free Path

Using the particle trajectory data sets from Table 1, it is possible to evaluate general properties that the ensemble of particles have after a sufficiently long duration of propagating in the turbulent magnetic field. The mean free path is the distance at which the instantaneous particle directions become uncorrelated with respect to those at time t=0. At this distance, and associated to scattering time scale, particles have lost memory of the direction distribution at initial conditions. This is a cumulative property of all particles, and it can be estimated by calculating the mean distance at which the direction of each particle has an angle of 90 from its initial position. This definition is equivalent to evaluating the velocity autocorrelation function and estimating the distance at which it is sufficiently close to zero (i.e., correlation to initial condition is lost).

The results of calculations are given in Table 1. As expected in this regime, the mean free path increases with energy. For the interarm region, where the injection scale is on the order of 100 pc, 60 pc for 30 PeV protons. For an energy four times smaller, 7.5 PeV, the mean free path decreases to a value of 50 pc. In the spiral arm, with injection scale on the order 10 pc, our calculation for 750 TeV protons is 5.0 pc. Note that in the two lowest energy data sets considered here, the particle gyroradius is on the order of the damping scale of the turbulent field. Therefore, pitch angle scattering arising from smaller scale magnetic perturbations is significantly reduced, and may result in the mean free path being overestimated to some degree. In the present work, the intent is to show the effects of turbulence in the particle spatial distribution in relation to the mean free path; therefore, whether the value is strictly correct is of limited importance in this context. In future studies, the need to use MHD simulations in a wider inertial range will be considered in more detail.

### 4.2. Sky Maps of Arrival Direction Distribution

The particle trajectories numerically integrated with Eqs. 1-2 for the sets in Table 1 are used to study the properties of their arrival direction distribution that results from the scattering off magnetic turbulence from a particle density gradient. As described in Section 2.2, the procedure used for determining the sky maps of the particles’ arrival distribution makes use of the trajectory integration backward in time, uniformly emanating from one point assumed to be Earth, until they exit the integration box. At a sufficiently long distance from the origin, particle trajectories accumulate in the direction of the mean magnetic field, since the perpendicular diffusion is smaller than that which is parallel to the magnetic field lines. A sphere of radius with center at the origin of the back-propagated trajectories is used to record the position and velocity direction of each trajectory. The trajectories are inverted in time from those locations on the sphere to the origin, by virtue of Liouville’s theorem (see Section 3 for a discussion of the validity of Liouville’s Theorem and the consequences of its use in the context of the present study).

It is likely that there are several sources of cosmic rays in the Milky Way, perhaps from different populations injecting particles into the ISM with different energy spectral shapes. It is also possible to suppose that a single source may dominate the cosmic ray distribution at Earth, depending on the energy range (see Section 1). In a scenario of isotropic diffusion, the cosmic ray density gradient is expected to be described by a simple dipole distribution. This is a similar distribution as would be expected if convective transport were the dominant source of the cosmic ray density gradient, such as in the case of superbubbles (Biermann et al., 2013) or the effect of Loop I shell in the local ISM (Schwadron et al., 2014). In the general and more likely scenario of anisotropic diffusion, cosmic ray propagation along the magnetic field line is faster than the perpendicular, thus producing a cosmic ray gradient ordered along the regular magnetic field (see, e.g., Kumar & Eichler (2014); Mertsch & Funk (2015)). Although the density gradient is not expected to be a dipole but rather a distribution that depends on the ratio of perpendicular to parallel diffusion, in this work it is assumed a simple dipole distribution, regardless of the origin of the underlying density gradient of the cosmic rays.

If the underlying distribution of cosmic rays is perfectly uniform, the effects of scattering off magnetic turbulence shuffles one isotropic distribution into another isotropic distribution. However, with an existing particle density gradient, the scattering processes redistribute particles from the region of higher density to that of lower density, and vice versa, thus creating a complex arrival distribution that can be described with higher order multipole terms of the spherical harmonic expansion.

In the process of inverting time on the numerically calculated trajectories, a dipole gradient distribution, assumed to be aligned with the direction of the mean magnetic field of the MHD snapshot, is introduced as a weight on each trajectory at the crossing point on a sphere of radius . The weight is calculated using the angle of the particle direction at the crossing point from that of the density gradient. The arrival distribution of these forward-propagated trajectories at Earth (i.e., the origin) is then determined. One key factor to remember is that only the small-scale angular anisotropy is studied, since this is the one that arises from the specific interaction with the turbulent magnetic field. Therefore, the large scale component, specifically the dipolar component of the map at Earth, is fitted and removed. Such a dipole component can have a different direction and amplitude than that injected. Figure 2 shows the sky map, in equatorial coordinates, of the arrival distribution of the 7.5 PeV protons at Earth before (on the left) and after (on the right) dipole subtraction, thus emphasizing the small scale features. A Gaussian smoothing with was used. The residual map shows medium- and small-scale angular structures arising from the breakout of the underlying dipole anisotropy after a propagation of = 50 pc.

The maps were created with the HealPix mapping tool (Gorski et al., 2005), which divides the sphere into pixels of equal areas. For the present work, a pixelization parameter of was used, which corresponds to 3072 pixels in total, with a mean spacing of . In Figure 2, the excess regions, with respect to the average isotropic flux, are identified by a red color and deficit regions by a blue one. Therefore, a pixel in which many particles pass through will be represented in a stronger red color than one that has only a few events.

Figure LABEL:fig:progression, under the assumption of the Earth’s location in the interarm zone with =100 pc, shows the sky maps progression of the 30 PeV cosmic ray arrival direction distribution with the dipole density gradient weight at different propagation distances of = 10 pc, 20 pc, 60 pc, and 90 pc. On each map, a dipole fit was performed, and the resulting dipole component was subtracted. Such a component may be different for each map, since it depends on the actual magnetic field structure at a given distance. The sky maps show that by increasing the propagation distance the arrival direction distribution progressively develops smaller structures up the mean free path (60 pc in this case, see Table 1). At larger distances, the arrival direction distribution reaches a statistically steady configuration (in Figure 2 only the 90 pc propagation distance is shown). Only propagation processes within the mean free path are responsible for the actual arrival direction distribution of cosmic rays at Earth, as discussed also in Giacinti & Sigl (2012). Whatever happens at larger distances is reshuffled by the scattering processes and is only relevant in the generation of the seed particle density gradient at large scale.

The actual distribution of cosmic rays depends on the specific realization of the turbulent magnetic field and on the particle energy. Figure 4 shows the steady-state arrival direction distribution of 750 TeV (on the left) and 30 PeV (on the right) cosmic rays. This figure qualitatively shows that in the higher energy case the anisotropy of the distribution shows more small-scale angular regions than in the lower energy case. For a quantification of such a visual property, it is necessary to calculate the angular power spectrum. The map from 7.5 PeV presents almost the same distribution as the one for 750 TeV, since the only differences between them are the assumption on the injection scale and that they are independent sets (as described in Section 2.2 and shown in Table 1). The particles at 7.5 PeV energy and are physically equivalent to particles at 750 TeV with .

### 4.3. Angular Power Spectrum

The sky maps of arrival direction distributions described in the previous sections are not dissimilar to experimental observations. However, it is the determination of the angular power spectrum that makes it possible to quantify the formation of small-scale structure anisotropy arising from scattering off magnetic turbulence. With the power spectrum, a sky map of arrival direction distribution is decomposed into spherical harmonics to provide information on the angular scale contribution of the anisotropy. The spectrum quantitatively indicates which multipole moments in the spherical harmonic expansion contribute to the observed map. The IceCube observatory provided a power spectrum of their Southern Hemisphere observation in the 10 TeV scale (Abbasi et al., 2011; Aartsen et al., 2015) and the HAWC gamma-ray observatory provided one for the Northern Hemisphere in the TeV scale (Abeysekara et al., 2014). The angular power spectrum is determined using the anafast tool in the HealPix framework.

Figure 5 shows the angular power spectrum from the numerical trajectory integration set of 30 PeV protons at propagation distances from the initial unperturbed dipole anisotropy at 10 pc, 20 pc, 60 pc and 90 pc, as in Figure LABEL:fig:progression, but without subtracting the dipole component. In the figure, the result from the observations made by the IceCube observatory is included as well. Note that the higher power values observed at low , not reproduced in the calculations, are most probably due to the partial sky view of the IceCube observatory, which causes correlations with the largest scale spherical harmonic moments. The experimental observation is for a median cosmic ray energy of about 20 TeV, much lower than the numerical calculations. The numerical calculation sets are normalized to the dipole component (i.e., = 1) of experimental power spectrum for the farthest propagation distance of 90 pc only. At shorter propagation distances, the normalization corresponds to the relative power obtained from the calculations. This normalization is valid since we are interested in the small-scale features, not on the assumptions on the large-scale anisotropy.

Note that the numerical calculation shows, within statistical fluctuations, that the power of the dipole component decreases with increasing propagation distance, due to the transfer of part of it into the higher components.
At very short distances, i.e., smaller than the mean free path, the low multipole moments are dominant. However, as the distance increases, more power is transferred to the higher multipole moments, which is the signature of particle interactions with the turbulent magnetic field. Once they reach the mean free path distance, in this case 60 pc, the power spectrum converges to a steady configuration, as the sky maps in Figure LABEL:fig:progression show.
In Figure 5, the power spectrum of an isotropic flux of the same number of particles is included as a gray band and properly normalized ^{1}^{1}1The power spectrum of the isotropic flux is calculated by generating 10,000 realizations of uniform particle direction distribution, matching the number of the integrated particle trajectories, and calculating the power spectrum for each of them. The 68%, 90%, and 99% containment bands are reported in Figure 5.
It is evident that the distribution at distance and beyond is not only a random distribution but also possesses a definite structure. It is noted that the angular power spectra calculated for propagation distance longer that 90 pc show the same relative normalization and shape of that corresponding to the mean free path. This is compatible with the achievement of steady state in the anisotropy structures as evident in the sky maps of Figure LABEL:fig:progression. The approach of the power spectra to the band of isotropic distribution for 30 means that, for the numerical realization studied here, the smaller angular scales are indistinguishable from random fluctuations of the isotropic flux.

As observed in Figure 4, the angular structure generated from the effect of scattering off magnetic turbulence depends on the particle energy (750 TeV and 30 PeV shown). This is evident also in the corresponding power angular spectra in Figure 6, calculated at the propagation distance of their mean free paths. The figure also shows the experimental result from the IceCube observatory at 20 TeV median energy, with the corresponding power of the isotropic distribution background and with the curve expected from the hierarchical breakup of angular components from Ahlers (2014), normalized to the dipole component of the IceCube result. Note that the angular power spectrum at higher energy is flatter than that at lower energy. This is compatible with the existence of more small-scale structures as evident in the sky maps of Figure 4.

## 5. Discussion

We have shown how small-scale anisotropy arises from the interaction of particles with the turbulent magnetic field. Specifically, we have addressed how the integration of trajectories in an MHD turbulent magnetic field provides a realistic understanding of the small-scale features present in the observations of cosmic ray anisotropy at Earth. The ISM is in a plasma state, where the MHD equations serve as a model for its dynamics, therefore an MHD turbulent magnetic field is a natural approach to study the magnetic field properties in the ISM. Previous work (Giacinti & Sigl, 2012) had considered the effects of magnetic turbulence on cosmic ray distribution. In this study, the authors have considered synthetic turbulence, which, on one hand, lacks the proper development of the gas dynamics but, on the other hand, provides the first qualitative connection between a magnetic turbulent field and cosmic ray arrival directions. In compressible MHD turbulence, scattering efficiency strongly depends on the wave type and how the particle gyroradius compares to the turbulence scales. Specifically, fast modes are identified as the main source for cosmic ray scattering (Yan & Lazarian, 2008).

The dynamics of the different turbulence modes and the relationship between particle energy and turbulence scale determine the actual scattering efficiency, which is most responsible for the breakout of a global cosmic ray density gradient into small-scale angular anisotropy. The angular power spectrum calculated for the data sets studied in this work appears to evolve in time until particles cross the mean free path distance. In a steady state, the shape of the angular power spectrum is a function of the magnetic field structure and of the consequent effects on particle propagation.

Studying particle trajectories in MHD magnetic turbulence provides a more realistic framework to investigate the behavior of cosmic ray propagation in the ISM, where turbulence is expected to be anisotropic, although MHD turbulence simulations typically represent a significantly more restrictive inertial range than the actual astrophysical plasmas.

The ISM is a complex environment and its exact representation is difficult to achieve; therefore, our MHD magnetic field can be considered one possible configuration of the magnetic field in the ISM. For this reason, direct comparison should be on the angular power spectrum itself, not on the exact topology of the small-scale features in our maps.

The framework for the application of Liouville’s theorem is provided in the context of cosmic ray arrival distribution and applied through the back-propagation method. Although particle trajectories appear to suffer from mild deviation from adiabaticity, Liouville’s theorem was applied in this particular case to study the first order effects of magnetic turbulence on the global topology of particle trajectories. This is because particles do not experience severe scattering in their trajectories, as shown in the first adiabatic invariant calculations; nonetheless, if the magnetic field were to vary dramatically with respect to the gyroradius of the particles, it prohibit application of Liouville’s theorem. The spread in the magnetic moment distribution in Figure 1 suggests that some trajectories may have experienced more scattering than average in their propagation. This effect will manifest in the anisotropy through higher power at high , since these particles will have had greater interactions with the turbulent field, resulting in a slightly flatter angular power spectrum. If the distribution on the plot had peaked at a higher value or a progressive trend had appeared on the mean magnetic moment plot of the data set, this would be an indication of a clear violation of Liouville’s theorem. Neither of these trends have statistically or significantly occurred in our calculations; nevertheless, we would expect them to appear in magnetic fields that have a strong variation in scale on the order of the gyroradius of the particles. Future studies will need to use MHD turbulence simulation with wider inertial range to include enlarge the energy range in which magnetic turbulence affects cosmic ray distribution. In this way, these studies will reveal the connection between the angular power spectrum of the cosmic ray arrival direction distribution and the turbulence properties, naturally accounting for spurious effects derived by the numerical methods used.

As mentioned in Section 2.2, the trajectory back-propagation method is intended to provide a high efficiency in the studies of particle propagation in magnetic fields. Such a numerical method provides a mapping of the regions in space that are magnetically connected to the arrival point, where particles are assumed to be isotropically distributed. Therefore, appealing to the conservation of the total power across all spherical harmonic contributions, as dictated by Liouville’s theorem, makes higher multipole moments possible. Anisotropy was studied as a function of an initial large-scale density gradient at a large distance from the arrival point. Such a density gradient, however, is generated by the same propagation processes that produce all angular structures in the arrival direction distribution as well. In fact, smaller scale structures can arise in trajectory forward-propagation integration methods, where particle escape is naturally accounted for, even without assuming an initial global anisotropy, as shown in Rettig & Pohl (2015). On the other hand, as long as a global anisotropy is developed at some distance larger than , independently on how it is generated, the arrival direction distribution reaches a steady state angular power spectrum, and the effects of anisotropy and observed anisotropy can be disentangled.

For the mean free path calculation, it is shown to be dependent on energy. In this regime, from 750 TeV to 30 PeV, the increases with energy. This is in agreement with how the more energetic particles interact with the perturbations of the magnetic field. In the case of the 750 TeV and 7.5 PeV sets, the minimum scale in the power spectrum of turbulence is of the size of the gyroradius of the particles, which may cause an overestimation of the , since there is limited power in the physical perturbations that interact with the particles. For our 30 PeV particles, we are well above this minimum scale, so the is unaffected.

The anisotropy maps and their corresponding angular power spectrum can tell us about the propagation of the particles themselves in the turbulent field. At very short distances, the low multipole moments are completely dominant and hold most of the power, as shown with the line of 10 pc in Figure 5. The reason for this is that the particles have not had enough time to interact with the features of the magnetic field, and the distribution is still reminiscent of the initial configuration. However, the particles continue to interact and structures start to develop. As we can see with the line of 20 pc in Figure 5, the higher multipole moments start to rise, with small-scale features becoming highly visible even in the maps in Figure 4. One interesting feature is that the small-scale structures develop within the mean free path, but once they reache this scale, the maps do not change significantly. This observation is confirmed in the angular power spectrum. Therefore, the distribution of power between the different multipole moments reaches a steady state. Of course, the particles keep moving and interacting after they have reached one , but the anisotropy itself does not change. From Figure 5 it is possible to see that this steady distribution is not isotropic, yet it possesses a definite structure that is dependent on the nature of the turbulent magnetic field. Consequently, the observed anisotropy is for the most part created in the last of the particle’s trajectory, and it becomes a way to indirectly probe the local ISM.

In Figure 5, we can see from the comparison with the observations that the experimental data behaves according to what we would have expected from a lower energy. In the case of 20 TeV, the distribution of power among the higher multipole moments is lower than at the higher energy, i.e., 30 PeV. One of the causes is that the 30 PeV particles interact with perturbations that carry more power than the ones at a lower energy. Therefore, when an interaction process occurs with these perturbations, the higher energy particles are affected more strongly and more small-scale structure is created.

This effect is confirmed in Figure 6, where the spectrum at the of energies 750 TeV and 30 PeV is compared with the observations at the median energy of 20 TeV. The distribution at 750 TeV is similar to the observations at 20 TeV; however, for 30 PeV, the distribution of high multipole moments is much higher and flatter than that at lower energy. Therefore, greater small-scale anisotropy is observed at higher energies. This is also easily seen in the maps in Figure 4.

The flatter angular power spectrum obtained with 30 PeV compared to 750 TeV protons tells us about differences in pitch angle scattering as a function of energy. However it is important to note that the 750 TeV data set corresponds to a gyroradius scale close to the dissipation region of the MHD turbulence inertial range. It is likely that scattering is underestimated, thus preventing a full development of small-scale angular structures in the particle anisotropy. From this point of view, it is reasonable to assume that the resemblance of the 750 TeV proton power spectrum to experimental data should be considered coincidental but still in agreement with the fact that lower energies should have less structure, as mentioned above. On the other hand, the difference of the 30 PeV proton power spectrum with experimental data does not necessarily mean that the fundamental processes responsible for the small-scale anisotropy are overestimated in the present study. The results show an energy dependence in the shape of the angular power spectrum that needs further study, requiring MHD turbulence simulation with a significantly wider inertial range.

Figure 6 includes the angular power spectrum resulting from hierarchical decay of angular scales if Liouville’s theorem is satisfied, as calculated by Ahlers (2014). Here the different curves have been normalized to the dipole of the observational data, as discussed in Section 4.3, so that only the small-scale structure is relevant. In Ahlers (2014), an existing global dipole anisotropy evolves to create higher multipole moments. The difference between that result and the one shown in the present work resides in the fact that here the shape of the distribution is determined by the specific turbulence characteristics of the magnetic field and the gyroradius of the particles. The 750 TeV case appears similar to Ahlers (2014), but since this set is at the damping scale and scattering is likely to be underestimated, we were already expecting less structure to be present.

Ahlers & Mertsch (2015) studied the effect of relative diffusion, i.e., from correlated nearby trajectories, as the contribution to the development of small scale anisotropy structures. In this complementary approach, the angular power spectrum is calculated based on the effect that a particle density gradient has on the trajectory topology shaped by homogeneous isotropic magnetic turbulence. In the present work, the relationship between the shape of angular power spectrum and the scattering processes induced by Alfvénic, slow, and fast modes in a compressible MHD turbulence were studied. This will make it possible to identify physical properties that shape the angular power spectrum, so that the problem can be inverted, i.e., the angular power spectrum and different cosmic ray energies and masses can be used to probe the properties of the interstellar magnetic turbulence.

## 6. Conclusions

This work explores the possibility that the local interstellar magnetic field could shape the arrival direction distribution of high-energy cosmic rays. An MHD turbulent magnetic field was used as a propagating medium that resembles the ISM, and particle trajectories with various energies were integrated in this field. To obtain the anisotropy in arrival direction distribution at Earth, the theoretical framework for application of Liouville’s theorem was provided, and the theorem was shown to be applicable in this specific case. Nonetheless, there could be cases where the magnetic field varies abruptly, which would ultimately violate the application of the theorem.

The results presented in this work show that small-scale anisotropy arises from the interaction of cosmic rays with the local turbulent magnetic field, as confirmed in the sky maps and angular power spectra. The angular power spectrum becomes flatter the higher the energy; therefore, experimental data at the 10 TeV scale is expected to be steeper than the numerical calculations presented in this paper. Cosmic rays with less rigidity are more sensitive to the lower power small-scale magnetic perturbations. Since cosmic rays are not composed only of protons, the contribution of heavy ions would yield a steeper angular power spectrum.

The inertial range of our turbulent magnetic field provides a limitation on the lowest energy that could be studied, in this case 750 TeV. Therefore, in the future, it will be necessary to extend this inertial range and sample even lower energies, so that direct comparisons with the observations at 20 TeV will be possible. Still, it is expected that at these TeV energies the effects of the heliosphere should be present as well, which would provide a primary topic for future study and publication. Another factor to improve is the number of particles that are propagated. If we were to have at least events, it would be possible to obtain a clearer view of the distribution of cosmic rays at Earth. Investigating how the properties of our local turbulent magnetic field might influence TeV-PeV cosmic ray arrival direction distribution will provide the basis for further exploring the observed anisotropy, and it will open the doors for a better understanding of our local interstellar medium.

## Appendix A Numerical Approach and Accuracy

Particle trajectories are calculated by integrating the 6D set of equations of motion, eqs. 1 and 2. As stated in section 2.2, the integration is performed numerically using the Bulirsch-Stoer method, which is considered one of the best known integration algorithms satisfying both high accuracy and efficiency (Press et al., 1986) and widely applied in the literature (e.g., Giacalone & Jokipii (1999) and Xu & Yan (2013)). The Bulirsch-Stoer integration algorithm, is a known method for numerical calculation of ordinary differential equation solutions, that combines the so-called Richardson extrapolation (to improve the rate of convergence of a sequence) and the modified midpoint method (which advances a vector of dependent variables y(x) from a point x to a point x + H by a sequence of n substeps each of size h). The result is that Bulirsch-Stoer algorithm provides high accuracy with relatively low computational effort. The accuracy of the algorithm is further controlled, during the numerical calculation, by monitoring the local truncation error estimated at each time step. If the relative error is larger than the relative tolerance level of 10, the step size is adaptively reduced in order to limit the error accumulation in both momentum and spatial coordinates, across the maximum integration time used in this work (corresponding to no more than 10,000 gyrations). Such error accumulation needs to be monitored for each specific problem in which the integration algorithm is used. In particular, the accuracy on the spatial and momentum coordinates were studied.

The performance of the Bulirsch-Stoer algorithm on the accuracy of momentum coordinates in our numerical calculation can be seen in Figure 7, where the particle energy variation (due to loss of accuracy) from that at is plotted as a function of the number of gyro-orbits (for a single particle and for the average of all particles used in our 30 PeV sample).

The maximum relative mean deviation from perfect energy conservation is found to be about 8.510 for the sample used in this work (100,000 particles in the 30 PeV energy set), although a single particle can reach a violation at the 1.610 level. This precision level in the energy conservation guarantees that particle trajectories are not significantly affected by numerical accuracy limitations, which are found to be marginal under the conditions of this study.

Since the adaptive time step algorithm constrains both spatial and momentum coordinates to the same relative error level, the accuracy in spatial coordinates is , even after 10,000 gyrations. Numerical diffusion, therefore, is limited to a level much smaller than Bohm diffusion, which is several orders of magnitude below diffusion induced by the stochastic wandering of magnetic field lines at all scales (Lazarian & Yan, 2014).

The effect on the particle set size was assessed in Xu & Yan (2013), which used the same integration stepping algorithms as in this work. In that paper, the perpendicular diffusion coefficient becomes stable when the sample size reaches about 1000 particles. The sets used in this study contain 100,000 or more particles (see Table 1), thus minimizing statistical accuracy effects on the global behavior of the ensemble of particles.

## References

- Aartsen et al. (2013) Aartsen, M. et al. 2013b, Astrophys. J. 765, 55
- Aartsen et al. (2015) Aartsen, M. et al. 2015, submitted to ApJ
- Abbasi et al. (2010) Abbasi, R. et al. 2010a, Astrophys. J. 718, L194
- Abbasi et al. (2011) Abbasi, R. et al. 2011a, Astrophys. J. 740 16
- Abbasi et al. (2012) Abbasi et al. 2012b, Astrophys. J. 746, 33
- Abdo et al. (2008) Abdo, A.A. et al. 2008, Phys. Rev. Lett. 101, 221 101
- Abdo et al. (2009) Abdo, A.A. et al. 2009, Astrophys. J. 698, 2121
- Abeysekara et al. (2014) Abeysekara, A.U. et al. 2014, Astrophys. J. 796, 108
- Aglietta et al. (2009) Aglietta, M. et al. 2009, Astrophys. J. 692, L130
- Ahlers (2014) Ahlers, M. 2014 Phys. Rev. Lett. 112, 021101
- Ahlers & Mertsch (2015) Ahlers, M., & Mertsch, P. 2015, arXiv:1506.05488
- Amenomori et al. (2005) Amenomori, M. et al. 2005, Astrophys. J. Lett. 626, L29
- Amenomori et al. (2006) Amenomori, M. et al. 2006, Science, 314, 439
- Amenomori et al. (2011) Amenomori, M. et al. 2011, Proc. 32nd ICRC, Beijing China
- Amenomori et al. (2007) Amenomori, M. et al. 2007, Proc. 30th ICRC, Mérida, Mexico
- Bartoli et al. (2013) Bartoli, B., et al. 2013, Phys. Rev. D 88-8, 082001
- Bartoli et al. (2015) Bartoli, B., et al. 2015, Astrophys. J. 809, 90
- Baumjohann et al. (1996) Baumjohann, W., & Treumann, R. A. 1996, Basic Space Plasma Physics (London: Imperial College Press)
- Beresnyak et al. (2011) Beresnyak, A., Yan, H., & Lazarian, A. 2011, Astrophys. J. 728, 60, 8 pp.
- Biermann et al. (2013) Biermann, P.L., Becker Tjus, J., Seo, E.-S., & Mandelartz, M. 2013, Astrophys. J. 768, 124
- Blasi & Amato (2012) Blasi, P., & Amato, E. 2012, JCAP 1, 11
- Bradt (2008a) Bradt, H. 2008, Astrophysics processes, Cambridge Univ. Press, Cambridge
- Bradt (2008b) Bradt, H. 2008, in Astrophysics Process (Cambridge, (supplement): Cambridge Univ. Press)
- Burkhart et al. (2014) Burkhart, B. et al. 2014, Astrophys. J. 790 130
- Cho & Lazarian (2002) Cho, J., & Lazarian, A. 2002, Phys. Rev. Lett. 88, 245001
- de Jong et al. (2011) de Jong, J. et al. 2011, Proc. 32nd ICRC, Beijing, China
- Desiati & Lazarian (2012) Desiati, P. & Lazarian, A. 2012, NPG 19, 351
- Desiati & Lazarian (2013) Desiati, P. & Lazarian, A. 2013, Astrophys. J. 762, 44
- Desiati & Zweibel (2014) Desiati, P. & Zweibel, E.G. 2014, Astrophys. J. 791, 51
- Eyink et al. (2011) Eyink, G.L., Lazarian, A. & Vishniac, E.T. 2011, Astrophys. J. 743, 51
- Effenberger et al. (2012) Effenberger, F. et al. 2012, A&A 547, A120
- Erlykin & Wolfendale (2006) Erlykin A.D., & Wolfendale A.W. 2006, Astropart. Phys. 25, 183
- Florinski et al. (2013) Florinski, V., Jokipii, J. R., Alouani-Bibi, F., le Roux, J. A. 2013, Astrophys. J. 776, L37
- Florinski et al. (2015) Florinski, V., Stone, E.C., Cummings, A.C., & le Roux, J.A. 2015, Astrophys. J. 803 47
- Frisch et al. (2012) Frisch, P.C. et al. 2012 Astrophys. J. 760, 106
- Frisch et al. (2014) Frisch, P.C. 2015 J. Phys.: Conf. Ser. 577, 012010
- Gaensler et al. (2011) B. M. Gaensler, et al. 2011, Nature 478, 214
- Gaisser et al. (2013) Gaisser, T.K. et al. 2013 Front. of Phys. 8-6, 748
- Giacinti & Sigl (2012) Giacinti, G., & Sigl, G., 2012 Phys. Rev. Lett. 109, 071101
- Giacalone & Jokipii (1999) Giacalone, J. & Jokipii, J.R. 1999, Astrophys. J. 520, 204
- Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, Astrophys. J. 438, 763
- Goldstein et al. (2002) Goldstein, H., Poole, C., & Safko, J. 2002, Classical Mechanics (3rd ed.; San Francisco, CA: Addison-Wesley)
- Gorski et al. (2005) Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
- Guillian et al. (2007) Guillian, G. et al. 2007, Phys. Rev. D 75, 062003
- Hall et al. (1999) Hall, D.L. et al. 1999, J. of Geophys. Res. 104, 6737
- Haverkorn et al. (2008) Haverkorn, B., et al. 2008, Astrophys. J. 680, 362
- Kumar & Eichler (2014) Kumar, R., & Eichler, D. 2014 Astrophys. J. 785, 129
- Jansson & Farrar (2012a) Jansson, R., & Farrar, G.R. 2012a, Astrophys. J. 757, 14
- Jansson & Farrar (2012b) Jansson, R., & Farrar, G.R. 2012b, Astrophys. J. 761, L11
- Jokipii (1966) Jokipii, J.R. 1966, Astrophys. J. 146, 480
- Jokipii & Parker (1969) Jokipii, J.R. & Parker, E.N. 1969, Astrophys. J. 155, 777
- Lazarian (2006) Lazarian A., 2006 AIPC, 874, 301
- Lazarian (2007) Lazarian, A. 2007, J. Quant. Spectros. & Radia. Transfer 106, 225
- Lazarian & Desiati (2010) Lazarian, A. & Desiati, P. 2010, Astrophys. J. 722, 188
- Lazarian & Vishniac (1999) Lazarian, A., & Vishniac, E.T. 1999, Astrophys. J. 517, 700
- Lazarian & Yan (2014) Lazarian, A., & Yan, H. 2016, Astrophys. J. 784, 38
- Lazarian et al. (2012) Lazarian, A. et al. 2012, Space Science Reviews 173, 557
- Lithwick & Goldreich (2001) Lithwick, Y., & Goldreich, P. 2001 Astrophys. J. 562, 279
- López-Barquero et al. (2016) López-Barquero, V. et al. 2016, to be submitted
- Mertsch & Funk (2015) Mertsch, P. & Funk, S. 2015 Phys. Rev. Lett. 114, 021101
- Minnie et al. (2009) Minnie, J., et al. 2009, JGR 114, A01102
- Munakata et al. (2010) Munakata, K. et al. 2010, Astrophys. J. 712, 1100
- Nagashima et al. (1988) Nagashima, et al. 1998, J. of Geophys. Res. 1031, 17429
- Manuel et al. (2014) Manuel, R., Ferreira, S., & Potgieter, M. 2014, Solar Physics 289, 2207
- Pogorelov et al. (2006) Pogorelov, N.V., Zank, G.P., & Ogino, T. 2006, Astrophys. J. 644, 1299
- Pohl & Eichler (2013) Pohl, M. & Eichler, D. 2013, Astrophys. J. 766, 9
- Ptuskin (2012) Ptuskin V. 2012, Astropart. Phys. 39, 44
- Press et al. (1986) Press, W.H., Flannery, B.P., & Teukolsky, S.A. 1986, Cambridge: University Press, 1986
- Rechester & Rosenbluth (1978) Rechester, A.B. & Rosenbluth, M.N. 1978, Phys Rev Lett 40, 38
- Rettig & Pohl (2015) Rettig, R. & Pohl, M. 2015, in Proc. of 34th ICRC, The Hague, The Netherland
- Santander et al. (2013) Santander, M. et al. 2013, in Proc. of 33rd ICRC, Rio de Janeiro, Brazil; arXiv:1309.7006
- Savchenko et al. (2015) Savchenko, V., Kachelrieß, M. & Semikoz, D.V. 2015, DOI:10.1088/2041-8205/809/2/L23
- Shuwang et al. (2011) Shuwang, C. et al. 2011, Proc. 32nd ICRC, Beijing China
- Shalchi (2009) Shalchi, A. 2009 Nonlinear Cosmic Ray Diffusion Theories, ASSL 362 (Springer)
- Schwadron et al. (2014) Schwadron, N.A., Adams, F.C., Christian, E.R., Desiati, P., Frisch, P., Funsten, H.O., Jokipii, J.R., McComas, D.J., Moebius, E., Zank, G.P. 2014, Science 343, 988
- Sridhar & Goldreich (1994) Sridhar, S., & Goldreich, P. 1994, Astrophys. J. 432, 612
- Sveshnikova et al. (2013) Sveshnikova, L.G. et al. 2013, Astropart. Phys. 50, 33
- Xu & Yan (2013) Xu, S., & Yan, H. 2013, Astrophys. J. 779, 140
- Yan & Lazarian (2008) Yan, H., & Lazarian, A. 2008, Astrophys. J. 673, 942
- Yan & Lazarian (2011) Yan, H., & Lazarian, A. 2011, Astrophys. J. 731, Issue 1, article id. 35, 10 pp. (2011)
- Zhang et al. (2009) Zhang, J.L. et al. 2009, Proc. 31st ICRC, Łódź, Poland
- Zhang et al. (2014) Zhang, M., Zuo, P., & Pogorelov, N. 2014, Astrophys. J. 790, 5