Small-scale Anisotropies of Cosmic Rays from Turbulent Flow

Within the classical convection–diffusion approximation, we show that the angular distribution of cosmic rays (CRs) in a highly turbulent flow may exhibit significant small-scale anisotropies. The CR intensity angular power spectrum C ℓ is then a direct reflection of interstellar turbulence, from which one expects C ℓ ∝ ℓ −γ−1 for ℓ ≫ 1, where γ is the power-law turbulence spectral index. Observations by IceCube and HAWC at TeV energies can be explained approximately with the Kolmogorov law γ = 5/3 with a convection velocity dispersion of 20 km s−1 on the scale of 10 pc.


Introduction
Cosmic rays (CRs) are mainly relativistic nuclei, with an energy spectral index (−∂ ln f /∂ ln p − 2) about 2.7 and relative intensity anisotropy on the order of 0.1% at TeV energies.The high level of isotropy implies that CRs are diffusive due to uniform pitch-angle scattering by magnetic field irregularities.This is consistent with the fact that the anisotropy is dominated by a dipole signal, which is believed to be associated with the spatial diffusion flux, and can thus be used to trace CR sources (Ahlers, 2016;Zhao et al., 2022;Zhang et al., 2022b;Qiao et al., 2023).
It is not surprising that the observed CR anisotropy deviates from a pure dipole (ℓ = 1), but the multipoles with spherical harmonic degrees ℓ > 1 cannot be explained by the standard diffusion theory.Besides some exotic scenarios (Kotera et al., 2013), most existing models for this issue focus on improving the statistical description of charged particles interacting with a classical electromagnetic field.It has been suggested, e.g., that the basic dipole configuration can be distorted by nonuniform pitch-angle scattering (Malkov et al., 2010) and/or non-diffusive transport (Harding et al., 2016), which may result from the local magnetic field structure (Giacinti & Sigl, 2012;Schwadron et al., 2014).
Interestingly, CR observations by the IceCube experiment and High-Altitude Water Cherenkov (HAWC) Observatory show that the TeV intensity sky map complies with a power-law angular power spectrum C ℓ ∼ C 1 ℓ −3 (Aartsen et al., 2016;Abeysekara et al., 2014Abeysekara et al., , 2018Abeysekara et al., , 2019)).This unique feature could be the key to understand the small-scale anisotropies, and has been interpreted as a consequence of the angular correlation of (two) particle trajectories with similar initial conditions, i.e., the particles undergo so-called relative diffusion in a turbulent magnetic field (Ahlers, 2014;Ahlers & Mertsch, 2015;Kuhlen et al., 2022).
In this paper, we shall propose an alternative point of view that ascribes C ℓ to turbulent convection under the classical diffusion approximation.Inspired by Earl et al. (1988), in previous work we have emphasized that nonuniform convection can induce dipole and quadrupole (ℓ = 2) anisotropies proportional to the diffusion coefficient, via inertial and shear forces, respectively.In particular, the shear anisotropy due to Galactic differential rotation may be important for PeV-EeV CRs (Zhang et al., 2022a).This paper aims to extend the convection scenario into the ℓ > 2 region.

Particle Back-Tracking
The basic assumption of the classical convection-diffusion approximation is that particles tend to completely lose information about their initial states due to scattering in the rest frame of the scattering center.In general, the scattering centers at different locations in a flow field have complicated relative motion, i.e., the convection velocity u is a function of the position vector r.For nonrelativistic convection, the particle momentum in the fluid rest frame is p = P − P u (r) /V , where P and V are the particle momentum and velocity in the observer's frame, respectively.Treating r and P as independent canonical variables, when a particle moves between neighboring locations separated by the mean free path λ, it obtains a momentum where ∆ denotes the difference between the neighboring systems, and v is the particle velocity in the fluid rest frame.Note that this expression excludes terms on O u 2 concerning inertial and relativistic corrections for convection.
For brevity, let us adopt the back-tracking argument (Ahlers & Mertsch, 2017;Zhang et al., 2022a), which embodies the physics of most interest.We need to work in the fluid rest frame, only in which can complete relaxation due to scattering be observed.Since Liouville's theorem requires statistical fluctuations to be balanced by kinematic terms, if the particle phase-space distribution approaches an isotropic state f after relaxation (so p∂f /∂p = p∂f /∂p), the unrelaxed fluctuation (assumed to be small) can be written as the time reversal increment ∆f = ∆ ICG f − λ • ∇f , where the second term corresponds to the diffuse dipole anisotropy, and the part associated with nonuniform convection, may be referred to as an inhomogeneous Compton-Getting (CG) effect.The conventional CG effect arises from uniform convection, which yields an additional dipole term ∆ CG f = (−P • u/V ) ∂f /∂P in the observer's frame.It is worth mentioning that the fluctuation continuously generated by a system under the fluctuation-relaxation equilibrium should be equivalent to ∆f .
As ∆u includes isotropic parts of the fluid deformation (e.g., the isotropic compression or expansion), we have ∆ ICG f = 0, where the angle brackets denote averaging over all directions.The fully random fluctuation arising from nonuniform convection should be ∆ ICG f − ∆ ICG f .See Appx.A for a more careful derivation with Bhatnagar-Gross-Krook (BGK) analysis.

The Role of Turbulence
Eq. (2) establishes a direct link between the microscopic anisotropy and macroscopic deformation.The small-scale fluctuations are terms associated with the tensor product of multiple directional vectors in the momentum space (λ/λ = p/p).Therefore, it is nontrivial to expand the system into a Taylor series around λ = 0, where U l−1 determines the lth term in the corresponding expansion of ∆ ICG f .On the other hand, ∆ ICG f can also be written as a spherical harmonic expansion, whose ℓth term defines the irreducible 2 ℓ -pole anisotropy.
The two expansions are generally inequivalent (see Appx.A; note the distinction between l and ℓ), while CR anisotropies are expressed via the spherical harmonics.
Here, the significance of the Taylor expansion is that it facilitates our understanding to the role of fluid regularity.In principle, the fluid description requires Λ λ such that the local system of particles can be relaxed, where Λ is the macroscopic scale on which the flow property changes significantly.Assuming u ∝ Λ ν , we make the following order-of-magnitude estimate, where Γ denotes the gamma function.This reduces to The factor (λ/Λ) l implies that the large-l corrections of ∆ ICG f are not important for Λ > λ.Considering that λ typically increases with the particle energy, in a regular nonuniform flow there should be no significant l > 2 anisotropy arising from convection, except for the most energetic particles with λ ∼ Λ.However, CR small-scale anisotropies are observed at different energies.To avoid the high-l cutoff over a broad spectrum, we need to consider a highly turbulent flow, which varies irregularly between any two spatial points so that Λ ∼ λ can be maintained at any energy.It can be seen from the detracing projection (Applequist, 1989) that the lth Taylor term is composed of 2 l−2n | 0 2n l, n ∈ Z -pole structures, hence we are led to conclude that turbulence is capable of producing polyenergetic multipole anisotropies that do not exponentially decay with ℓ.
One may ask: Since structures with Λ < λ are hydrodynamically unresolvable, how can such small-scale anisotropies have a fluid description?To answer this, it is necessary to distinguish between fore-and back-ground flows: the former represents the flow of particles under consideration, while the latter is that of scattering centers.The establishment of relaxation means that the two flows share the same u on the scale λ, whereas the background flow is also defined on smaller scales.It is ∆u, which can now be interpreted as the background nonuniformity, and thus as a superposition of the background small-scale structures, that determines the configuration of the foreground small-scale anisotropies.
The spectrum of the structures is encoded in the Fourier expansion of the velocity correlation function.But we are more interested in the multipole expansion where Y LM is the spherical harmonic function.Using the plane-wave expansion (Landau & Lifshitz, 1991), one has 3 , where j L refers to the spherical Bessel function of the first kind, and ũ is the Fourier transform of ∆u.
Let the overline represent an ensemble average over all possible realizations of the local field.As the two-point correlation ∆u (λ (Landau & Lifshitz, 1987).For convenience, let us consider an equipartition system with ũα , where δ αα ′ is the Kronecker delta (α and α ′ = x, y or z), and w is the omnidirectional turbulence spectrum.Then, from the orthogonality condition 4π (Kuhlen et al., 2022), we have and implying that the turbulence energy is shared equally among all directional degrees of freedom.Note the identity Feldbrugge, 2023, who used χ L = 1/2), and we may take the approximation In this regard, there is a simple geometric interpretation: Structures described by the Lth term of Eq. ( 5) have an effective angular scale 2π/ (L + χ L ), corresponding to two points separated by a length scale 2πλ/ (L + χ L ) on a sphere of radius λ, while in the interval k L+1 − k L the velocity variance between the two points is roughly u 2 L .So the dependence of u 2 L on L is similar to that of w (k L ) on k L ; whereas they are different by a smoothing factor since the width of j 2 L is actually nonzero.As an example, the exact result of Eq. ( 6) for a power-law spectrum w ∝ k −γ with γ < 2L + 1 is ] in the large-L limit.

Multipole Anisotropies
According to Eqs. ( 2) and ( 5), the 2 ℓ -pole moment of ∆ ICG f is determined by the dipole-2 L -pole coupling where b LM,±1 = 2π/3 (ic LMy ∓ c LMx ), b LM0 = 4π/3c LMz , and Selection rules of the Wigner 3j symbols (Landau & Lifshitz, 1991) 6), ( 10) and ( 11).The cyan dot-dashed curve marks the Ahlers (2014; A14) model.The green short-dashed curve refers to modeling by Ahlers & Mertsch (2015;AM15).The black solid line shows the noise level dominated by HAWC, with the dark and light gray bands indicating the 68% and 95% confidence levels, respectively.Though the noise level for IceCube alone is not displayed in Abeysekara et al. (2019), from the sky coverage and number of events one can estimate it to be around 10 −11 (Ahlers & Mertsch, 2017).Right: Turbulence spectra used in the calculation.Note that here we set jL (kLλ) = max jL for labeling kL.
) δ mm ′ , while the angular power is determined by the diagonal elements which are m-independent.Therefore, any dependence of |a ℓm | 2 on m would imply that the equipartition theorem is broken in an asymmetric system.

Data Fit
It is mathematical and physical simplicity that motivates us to consider the equipartition hypothesis.Unfortunately, moduli of spherical harmonic coefficients of the TeV CR intensity have been reported to be mdependent (see Tab. 3 in Abeysekara et al., 2019).Nevertheless, we notice that the coefficient of variance for such moduli at any fixed ℓ is smaller than 1.Moreover, it remains possible that the observing time is insufficient.On the other hand, as Eq. ( 11) washes out the dependence on m, there must be an effective equipartition flow that produces the same C ℓ .For simplicity, let us assume the effective flow to be an acceptable approximation to the realistic one; especially, they have similar direction-averaged spectra.
In Fig. 1, we show that the observed CR angular power spectrum at a median energy of 10 TeV can be described for ℓ > 1 with the above turbulent convection scenario, which fulfills the Kolmogorov law γ = 5/3.The ℓ = 1 term is not included in our calculation because it involves the regular field u 0 (see Eq. 10).Given the higher noise level of HAWC, we fit only the IceCube data to cover the highest multipole moment provided (ℓ = 30), while the two experiments indicate basically the same spectral shape for 1 < ℓ 15.Assuming ∂ ln f /∂ ln p = −4.7, the data fit with Eq. ( 7) (i.e., Λ c → ∞) gives σ (10 TeV) ≈ 23 km/s, where σ = 2πw (2π/λ) /λ may be recognized as the CR convection velocity dispersion.For comparison, we also display the fitting results from other two standard scenarios concerning relative diffusion (Ahlers, 2014;Ahlers & Mertsch, 2015), whose analytic expressions can be found in Appx.B.
In fact, scattering centers of CRs are more likely to be plasma waves, whose propagation may be characterized with the Alfvén velocity v A = B/ √ 4πρ, where B is the magnetic field, and ρ is the mass density.Propagation of uncertainty yields ( 2 , with σ vA , σ B and σ ρ the corresponding rootmean-square values.On the other hand, the quasilinear diffusion theory claims where r g is the gyroradius, and σ 2 B (2πr g ) = σ 2 B (λ) (2πr g /λ) γg−1 is the mean-square magnetic field at the resonance wavenumber 1/r g with γ g the power-law spectral index of turbulence on the same scale (Strong et al., 2007).Setting λ = 10 pc, r g = 0.003 pc and γ g = γ = 5/3, for the local ISM with B ∼ 3 µG and ρ ∼ 0.5σ ρ (10 pc) ∼ 1.67×10 −25 g/cm 3 (Lee & Lee, 2019) we obtain σ vA (10 pc) ∼ 20 km/s, which is comparable with the fitted σ (10 TeV).
Consequently, we expect C ℓ (10 TeV) ∝ ℓ −8/3 for ℓ ≫ 1 until k L (10 TeV) reaches a broken wavenumber, around which the turbulence spectral index changes.Indeed, observational and theoretical uncertainties still admit a moderate variation of the index, while any scale dependence of the index must be imprinted on C ℓ .For example, as it has been claimed that the turbulence power may be enhanced (with respect to a single power-law spectrum) below the heliospheric scale of 2πη × 120 AU (η ∼ 1; Zank et al., 2019;Lee & Lee, 2020), we accordingly predict an enhancement of C ℓ for ℓ λ/ (120η AU), which is however far beyond the resolution of existing TeV observations.
With the same γ and σ (10 TeV) mentioned before, we plot energy variations of C ℓ in Figs.2-4.As seen, although turbulent convection with a power-law spectrum can provide a rough description to the angular power  spectra measured by HAWC (above the noise level) at multi-TeV energies, there is a tendency that the observed energy dependence in 10-100 TeV deviates from the power law below 10 TeV.Can this "anomaly" be ascribed to the scale dependence of the turbulence spectral index?Intuitively, there are two scenarios for the amplitude decrease with energy in 10-100 TeV: (i) γ > 1 and γ g > 2, which may partly account for efficient damping of plasma waves on scales between 2πr g (10 TeV) ∼ 0.02 pc and 2πr g (100 TeV) ∼ 0.2 pc; (ii) γ < 1 and γ g < 2, which may partly arise from energy injection by supernova remnants (SNRs) at a correlation length Λ c ∼ 50 pc.
In the following, we shall quantify only the injection effect from SNRs, and leave other effects for future work.For simplicity, let us consider the steady-state solution of the wave diffusion equation with monochromatic injection (Zhou & Matthaeus, 1990) , where Θ denotes the unit step function.The results in Fig. 4 suggest that the anomaly cannot be explained with the SNR injection effect alone, because, most importantly, the data exhibit no strong ℓ-dependence of the critical momentum p c , around which the energy dependence of C ℓ changes, while the theory predicts p c ∝ ℓ 1/(2−γg) according to k L (p c ) ∼ 2π/Λ c .Nevertheless, spectral flattening of SNR-driven turbulence on tens-of-pc scales is of physical significance (Chepurnov & Lazarian, 2010).In addition, the ℓ-independent p c may reflect a change in the energy dependence of λ, which is controlled by w (1/r g ).
Above 100 TeV, existing observations indicate that the energy dependence of the anisotropy amplitude projected onto the equatorial plane becomes a "standard" power law again, but the normalization factor is significantly smaller than that below 10 TeV (Zhang et al., 2022a).Recently, similar results of the angular power have also been reported preliminarily by IceCube using 11 years of data (McNally et al., 2023).We accordingly speculate that interstellar turbulence has a low-level component (compared with the SNR-driven one), which is injected due to Galactic large-scale convection at a correlation length no less than kpc.

Summary
With the idea that particle distributional fluctuations can be related to fluid local nonuniformity, we have demonstrated that turbulent convection may give rise to remarkable small-scale anisotropies of the distribution.Of particular note is that this scenario is based on the classical convection-diffusion approximation, which is ruled out by other models for the CR small-scale anisotropy problem.
Our model, described by Eqs. ( 6), ( 10) and ( 11), indicates that the 2 ℓ -pole anisotropy total power (2ℓ + 1) C ℓ is roughly proportional to w (ℓ/λ) /λ, the variance of the turbulent velocity at the scale 2πλ/ℓ in the characteristic wavenumber interval 1/λ.A power-law spectrum w ∝ k −γ exactly leads to C ℓ ∝ ℓ −γ−1 for ℓ ≫ 1.The fit to observational data brings an attractive picture that interstellar turbulence, with the Kolmogorov law γ = 5/3 valid at least on scales of 2-30 pc, and a velocity dispersion about 20 km/s for Alfvénic convection on the scale of 10 pc, is responsible for CR multipole anisotropies in the TeV range.This hypothesis remains to be further tested by future precise measurements of broad-spectrum anisotropies and turbulence.
where q is defined via the relative scattering rate ν r ∝ (1 − cos θ) q , with θ the angle between directions of motion of two particles.The fitting result is q ≈ 0.6, implying C ℓ ∝ ℓ −2.8 for ℓ ≫ 1 in AM15.
As C 1 is believed to be dominated by the diffuse dipole anisotropy, the energy dependence of C ℓ in A14 and AM15 can be affected largely by the distribution of CR sources, which is not the subject of this study.

Figure 1 :
Figure1: Left: Fits to the CR angular power spectrum at 10 TeV, with the observational data fromAbeysekara et al.  (2019).The red solid and blue long-dashed curves represent scenarios described by Eqs.(6), (10) and (11).The cyan dot-dashed curve marks the Ahlers (2014; A14) model.The green short-dashed curve refers to modeling byAhlers & Mertsch (2015; AM15).The black solid line shows the noise level dominated by HAWC, with the dark and light gray bands indicating the 68% and 95% confidence levels, respectively.Though the noise level for IceCube alone is not displayed inAbeysekara et al. (2019), from the sky coverage and number of events one can estimate it to be around 10 −11(Ahlers & Mertsch, 2017).Right: Turbulence spectra used in the calculation.Note that here we set jL (kLλ) = max jL for labeling kL.

Figure 2 :
Figure 2: Multi-energy C ℓ corresponding to fitting results in Fig. 1, with the observational data from Abeysekara et al. (2018; including systematic errors).The colored bands show the noises at a 90% confidence level.