Numerical Testing of Mirror Diffusion of Cosmic Rays

The tension between recent observations and theories on cosmic-ray (CR) diffusion necessitates exploration of new CR diffusion mechanisms. We perform the first numerical study on the mirror diffusion of CRs that is recently proposed by Lazarian & Xu. We demonstrate that the perpendicular superdiffusion of turbulent magnetic fields and magnetic mirroring that naturally arise in magnetohydrodynamic (MHD) turbulence are the two essential physical ingredients for the mirror diffusion to happen. In supersonic, subsonic, and incompressible MHD turbulence, with the pitch angles of CRs repeatedly crossing 90° due to the mirror reflection, we find that the mirror diffusion strongly enhances the confinement of CRs, and their pitch-angle-dependent parallel mean free path can be much smaller than the injection scale of turbulence. With the stochastic change of pitch angles due to gyroresonant scattering, CRs stochastically undergo slow mirror diffusion at relatively large pitch angles and fast scattering diffusion at smaller pitch angles, resulting in a Lévy-flight-like propagation.

In the direction parallel to the magnetic field, most earlier studies are focused on the diffusion induced by gyroreschao.zhang@ufl.edu;xusiyao@ufl.eduonant scattering.Due to the scale-dependent anisotropy of Alfvén and slow modes in MHD turbulence (Goldreich & Sridhar (1995); Lazarian & Vishniac (1999)), Alfvén and slow modes are inefficient in scattering the CRs with gyroradii much smaller than the turbulence injection scale (Chandran (2000); Yan & Lazarian (2002); Xu & Lazarian (2020)), while fast modes are identified as the more efficient agent of scattering (Yan & Lazarian (2004)).More recently, strong scattering of CRs by sharp intermittent magnetic field bends in MHD turbulence is proposed to be important for affecting CR parallel diffusion (Lemoine (2023); Kempski et al. (2023); Butsky et al. (2023)).The enhanced local scattering can be associated with the intermittent and fractal structure of MHD turbulence (Isliker & Vlahos (2003)).
The resonant scattering faces its long-standing 90 • problem in the framework of the quasi-linear theory (Jokipii (1966)).To resolve this issue, nonresonant interactions such as magnetic mirroring was explored (Fermi (1949); Cesarsky & Kulsrud (1973), CK73 henceforth;Noerdlinger (1968); Klepach & Ptuskin (1995); Chandran (2000)).This consideration of the mirroring effect naturally solves the 90 • problem and limits the pitch angle range for gyroresonant scattering (CK73; Xu & Lazarian (2020)).Magnetic mirrors can naturally form in MHD turbulence due to compressions of magnetic fields, which are induced by pseudo-Alfvénic modes in an incompressible medium and slow and fast modes in a compressible medium.Unlike the trapping effect of magnetic mirrors considered for compressible MHD waves in CK73, Lazarian & Xu (2021) (LX21 henceforth) found that in MHD turbulence due to the perpendicular superdiffusion of turbulent magnetic fields (Lazarian & Yan (2014)), CRs do not experience trapping.Instead, they bounce among different magnetic mirrors and move diffusively along the local magnetic field, which leads to a new diffusion mechanism termed "mirror diffusion".The mirror diffusion accounts for the suppressed diffusion of CRs when the mirroring condition is satisfied (Xu (2021); Barreto-Mota et al. (2023)).
In this work, we will carry out the first numerical test of the mirror diffusion with test particle simulations in MHD turbulence.We will also examine the pitch-angle-dependent mean free path of mirror diffusion that is analytically predicted by LX21.The outline of this letter is as follows.In Section 2, we describe the numerical methods used for the test particle simulations.In Section 3, we present our numerical results on testing the mirror diffusion and measuring its mean free path in different MHD turbulence regimes.Our conclusions can be found in Section 4.
2. NUMERICAL METHODS For the turbulent magnetic fields used for our test particle simulations, we take snapshots from 3D compressible and incompressible MHD turbulence simulations after the MHD turbulence reaches a statistically steady state.The compressible MHD turbulence simulations are taken from Hu et al. (2022), which are generated via ZEUS-MP/3D code (Hayes et al. (2006)).It solves the ideal MHD equations in a periodic box under the isothermal condition.The turbulence is solenoidally driven in Fourier space with the injected turbulent kinetic energy peaked at wavenumber k = 2, which is half of the data cube size, approximately 400 grids for the subsonic data cube (M1 henceforth) and 250 grids for the supersonic data cube (M2 henceforth).The numerical dissipation occurs at l diss ≈ 20 grids for M1 with numerical resolution = 7923 and l diss ≈ 25 grids for M2 with numerical resolution = 512 3 .The incompressible MHD turbulence simulation (M3 henceforth) is taken from Cho (2010), which is generated by using the pseudo-spectral code developed by Cho & Vishniac (2000).It solves the incompressible MHD equations in a periodic box with hyperviscosity and hyperdiffusivity.The peak of turbulent energy injection is at k ≈ 2.5 (corresponding to the turbulence injection scale L inj ≈ 200 grids), and the numerical resolution is 512 3 .The numerical dissipation occurs around l diss ≈ 10 grids.The inertial range of MHD turbulence is within [L inj , l diss ].In M1 and M2, the mean magnetic field ⃗ B 0 is in the z direction and in M3 ⃗ B 0 is in the x direction.By applying the hyperviscosity and hyperdiffusivity, M3 has a slightly more extended inertial range compared to M2 with the same resolution.Further details about these simulations can be found in Hu et al. (2022), Cho & Vishniac (2000) and Cho (2010).Since we are interested in the diffusion of relativistic CRs in non-relativistic MHD turbulence, that is, the test particle speed is much greater than the turbulent speed and the Alfvén speed V A , the turbulent magnetic fields can be treated as a static background for the test particle and the motional electric field can be neglected1 .The parameters of the three MHD simulations used in this work are summarized in Table 1.MHD turbulence is characterized by the Alfvén Mach number M A = V L /V A ≈ ⟨v⟩/ ⟨B⟩/ 4π⟨ρ⟩ , where V L is the turbulent velocity at L inj , ⟨v⟩, ⟨ρ⟩ and ⟨B⟩ are the rms values of velocities, densities, and magnetic field strength. 2The three MHD simulations we use are sub-Alfvénic with M A < 1. 3 For compressible MHD turbulence, M s = V L /c s is the sonic Mach number, where c s is the isothermal sound speed.We include both subsonic (M s < 1) and supersonic (M s > 1) cases for the compressible MHD turbulence.
For our test particle simulations, we inject relativistic test particles that represent CR particles in a snapshot of MHD turbulence simulations.To trace the test particle trajectory, we numerically solve the following equations: for a particle of position ⃗ r, velocity ⃗ u, charge q and mass m. ⃗ B is the magnetic field, c is the speed of light, and t is the time.We follow a numerical method similar to that adopted in Xu & Yan (2013); Hu et al. (2022).We use the Bulirsch-Stoer method to integrate Eq. 1 along with 2 with an adaptive time step (Press et al. (1988)).It provides a high accuracy for the particle trajectory integration with a relatively low computational effort (see also López-Barquero et al. (2016)).At each time step, the local magnetic field at the location of a test particle is interpolated by the cubic spline using the magnetic fields of the 10 3 neighbouring points.
The energy of a test particle is well conserved during the simulation and is quantified by the ratio r L /L inj , where the Larmor radius r L is defined as with the Lorentz factor γ. Throughout this paper, we use L inj to normalize length scales as L inj is relatively well determined in our simulations.The gyrofrequency is defined as Ω = u/r L .At every time step, we measure µ = cos θ, where θ is the angle between ⃗ u(⃗ r) and ⃗ B(⃗ r) , i.e., the pitch angle of a test particle.To achieve statistically reliable results, we inject a sufficiently large number of test particles and take an ensemble average of the numerical measurements over all test particles.

NUMERICAL RESULTS
The mirror diffusion of CRs in MHD turbulence arises from both the mirror reflection by magnetic compressions and the perpendicular superdiffusion caused by the Alfvénic turbulence.The latter plays an essential role in enabling the parallel diffusion of mirroring particles (LX21), which would otherwise be trapped between two fixed magnetic mirror points (CK73).To numerically test the mirror diffusion, we will first demonstrate the perpendicular superdiffusion of CRs, which has been analytically studied (Lazarian & Yan (2014); Lazarian et al. ( 2023)).The previous numerical testing of perpendicular superdiffusion of CRs was performed at small θ (e.g., Xu & Yan (2013); Hu et al. (2022)).To demonstrate the perpendicular superdiffusion of mirroring CRs with large θ (see Section 3.2.1),we will adopt large initial θ.

Perpendicular Superdiffusion of mirroring CRs
The perpendicular superdiffusion of CRs is caused by that of turbulent magnetic fields and takes place within the inertial range of MHD turbulence (Lazarian & Yan (2014)).To have initially small separations between the particle trajectories, we inject 50 beams of test particles, and each beam consists of 20 particles.Within each beam, the 20 particles are initially randomly distributed within a volume of size of around 5 to 8 grids, so that the average initial separation between a pair of particles is around 3 to 5 grids.All test particles have the same initial pitch angle cosine µ 0 ≈ 0.10 and r L = 0.01L inj .Within each beam, we measure the perpendicular separation δ x of every pair of particle trajectories with respect to ⃗ B 0 and obtain its rms value for all pairs.We then take the average over all beams as the measured perpendicualr displacement ⟨δ x2 ⟩ .and M3 as a function of t (normalized by Ω −1 ) is presented in Fig. 1.In all three cases, irrespective of the MHD turbulence regime, perpendicular superdiffusion is observed above l diss (see values in Table 1).Below l diss , CRs undergo subdiffusion with ⟨δx 2 ⟩ ∝ t α and α < 0.5.Above l diss , we see the transition from perpendicular superdiffusion with α > 0.5 to perpendicular normal diffusion with α ≈ 0.5.For the perpendicular superdiffusion, α = 0.75 is expected for the case when CRs undergo diffusive motion along the magnetic field (Lazarian et al. (2023)).Our result indicates that the mirroring CRs simultaneously undergo parallel diffusion and perpendicular superdiffusion within the inertial range of MHD turbulence.

Mirror Diffusion
In the presence of a magnetic mirror with a longitudinal magnetic gradient, CRs with r L smaller than the mirror size l ∥ and µ satisfying can be reflected by the magnetic mirror under the mirror force, where δB ∥ is the parallel magnetic fluctuation over l ∥ and B 0 is the mean magnetic field strength.The mirror force is (CK73) where M is the magnetic moment where u ⊥ is the perpendicular component of ⃗ u.For mirroring CRs with their motion along the magnetic field dominated by magnetic mirroring, there is no stochastic change of µ, and M is a constant, known as the first adiabatic invariant.As an example, in Fig. 2 we present the trajectory of one test particle with r L = 0.01L inj and µ 0 = 0.05 in M2.We adopt a sufficiently small r L (< L diss ) for mirroring to dominate over gyroresonant scattering and a sufficiently small µ 0 for the condition in Eq. ( 4) to be satisfied.The smallest r L is limited by our numerical resolution.In fact, r L = 0.01L inj is close to the smallest numerically resolvable scale, i.e., the grid size.In Fig. 2 (a), the particle trajectory is color coded by B/B 0 , where B is the magnetic field strength.The corresponding F ∥ points toward a region with a weaker B/B 0 .
It can be seen that the particle is reflected at different mirror points (i.e., regions with stronger B) multiple times and move back and forth in the direction parallel to the magnetic field.However, in the direction perpendicular to the ⃗ B 0 in the z direction, the particle is not trapped due to the perpendicular superdiffusion of turbulent magnetic fields (see Section 3.1).As expected, CRs are not trapped between two magnetic mirorr points, but diffusively move along the magnetic field in MHD turbulence due to the perpendicular superdiffusion.

Transition between mirror diffusion and scattering diffusion
CRs with sufficiently small µ (see Eq. 4) are subject to the nonresonant mirroring interaction and undergo the mirror diffusion.For CRs with larger µ, their parallel diffusion is dominated by the resonant pitch-angle scattering and undergo the scattering diffusion (LX21).The transition between mirror diffusion and scattering diffusion takes place as µ stochastically changes due to the gyroresonant scattering.According to LX21, the transition between the two diffusion regimes is expected to occur at a critical µ as µ c .In compressible MHD turbulence, under the consideration that fast modes are mainly responsible for mirroring and scattering, µ c is given by (LX21) where and corresponds to the pitch angle cosine where mirroring and scattering are in balance, with δB f as the magnetic fluctuation at L inj of fast modes.Via mode decomposition (see Cho & Lazarian (2003); Lazarian & Pogosyan ( 2012)), we find that δB f /B 0 is close to 0.40 for both M1 and M2, yielding µ max ≈ 0.52.For incompressible MHD turbulence, pseudo-Aflvén modes are responsible for mirroring.Under the consideration of inefficient scattering by Alfvén and pseudo-Alfvén modes for CRs with r L ≪ L inj , we have (LX21) where δB s is the magnetic fluctuation of pseudo-Alfvén modes at L inj .As pseudo-Alfven modes are slaved to Alfven modes (Beresnyak & Lazarian (2019) Next we numerically demonstrate the transition between the two different diffusion regimes.Here we choose a slightly larger r L compared to that in Fig. 2 for particles to undergo more scattering.As shown in Fig. 3, from the trajectory of a test particle with r L = 0.02L inj and µ 0 = 0.05 in M2, we can clearly see the transition from the slow mirror diffusion within a small region to the fast scattering diffusion over a large distance.During the mirror diffusion, the particle is confined in the direction of ⃗ B 0 (the z direction) with multiple mirror reflections.While during the scattering diffusion, the particle does not reverse its moving direction and rapidly moves along the magnetic field.
Furthermore, in Fig. 4 we present the time evolution of µ, particle position in the z direction, and 2M/mu 2 corresponding to the particle trajectory shown in Fig. 3.The numerical measurements and their gyro-averaged results are shown by blue and orange lines, respectively.Before Ωt ≈ 200, we see oscillations of µ around µ = 0, i.e., θ = 90 • , which corresponds to the mirror reflection at the mirror points.Accordingly, the particle moves back and forth in the z direction (middle panel in Fig. 4).After Ωt ≈ 200, µ changes stochastically, and the particle is not confined in the z direction.The transition between mirror diffusion and scattering diffusion occurs at µ ≈ 0.5, which is consistent with our estimate for µ max in M2 (see Eq. 8).As expected, M is conserved during the mirror diffusion, but varies dramatically during the scattering diffusion as efficient scattering violates the first adiabatic invariant (lower panel in Fig. 4).The distinctive features of mirror diffusion and scattering diffusion can also been seen in simulations with a lower resolution (see Appendix).
The transition between mirror diffusion and scattering diffusion happens due to the stochastic change of µ caused by scattering.To illustrate the stochastic change of µ, in Fig. 5, we present the time evolution of the probability density distribution of µ of 5000 particles with r L = 0.01L inj and µ 0 = 0.05 in M2.We see that due to the stochastic scattering, the distribution of µ gradually evolves into a more uniform distribution.As a result, CRs stochastically undergo the slow mirror diffusion when µ < µ c and the fast scattering diffusion when µ > µ c , leading to a Lévy-flight-like propagation.

Pitch-angle-dependent Parallel Mean Free Path of Mirror Diffusion
For the scattering diffusion, the scattering mean free path corresponds to the change of µ over the range [µ c , 1] (CK73; LX21).In the presence of the mirror diffusion at µ < µ c , the 90 • problem of scattering, i.e., the vanishing scattering near 90 • (Jokipii (1966)), is naturally solved, and the scattering mean free path becomes finite.For the mirror diffusion, its parallel mean free path λ ∥ is determined by the size of the magnetic mirrors that the mirroring particles mainly interact with (LX21).For CRs with µ < µ c , they preferentially interact with the mirrors satisfying (CK73; Xu & Lazarian (2020)) where l ∥ is the mirror size, and δb(l ∥ ) is the parallel magnetic fluctuation at l ∥ .Along the energy cascade of MHD turbulence, toward smaller scales, δb(l ∥ ) decreases, but the magnetic field gradient δb(l ∥ )/l ∥ increases, and the corresponding mirroring becomes more efficient (see Eq. 5).The above equation suggests that among the mirrors with sufficiently large magnetic fluctuations to reflect the CRs with µ, the smallest ones dominant the mirroring.Therefore, λ ∥ of mirror diffusion is expected to depend on µ (LX21).λ ∥ quantifies the slow diffusion of CR particles with sufficiently large pitch angles in the vicinity of CR sources (e.g., Xu (2021)).
To test the µ-dependence of λ ∥ , we inject test particles with different µ 0 .For the particles that exhibit oscillations within the same range [−µ, µ], we measure their displacement ∆z in the z direction between two mirror points (see Fig. 4 as an illustration) and take the average as the measured mean free path parallel to ⃗ B 0 corresponding to µ.For sub-Alfvénic turbulence, magnetic field lines are mildly perturbed by turbulence especially at small scales.Our measurement serves as a good approximation of ∥ parallel to the local magnetic field at small µ.
In Fig. 6, we present the measured λ ∥ as a function of µ for M1, M2, and M3 up to the theoretically estimated µ max (Eqs.( 8) and ( 10)).The width of each µ bin is 0.05, and the number of test particles in each µ bin ranges from 50 to 100 .We adopt r L = 0.01L inj ≲ l diss (see Table 1) for all the test particles to suppress the gyroresonant scattering, but the nonresonant mirroring is still active.Therefore, we have µ c ≈ µ max for both compressible and incompressible MHD turbulence.In all three MHD turbulence regimes, we see that λ ∥ increases with increasing µ within the range µ < µ max , as predicted by LX21.The error bar represents the standard deviation σ and also increases with increasing µ.At a given µ, we find that ∆z has a distribution (see Fig. 7).As the degree of field line wandering increases toward larger scales, the deviation of the local magnetic field direction from the mean magnetic field direction becomes more significant at large scales.This may give rise to a larger σ of the ∆z distribution at a larger µ (see Fig. 7).
4. CONCLUSIONS With test particle simulations in MHD turbulence, we performed the first numerical test of mirror diffusion of CRs.
Here we summarize our main findings.
(1) We identify the CRs with µ < µ c as mirroring CRs, which are subject to the mirror reflection by magnetic compressions in both compressible and incompressible MHD turbulence.We find that the mirroring CRs simultaneously un- field and mirror diffusion along the turbulent magnetic field.Rather than being trapped between two fixed magnetic mirror points, the perpendicular superdiffusion enables encounters of mirroring CRs with different magnetic mirrors and their diffusion parallel to the magnetic field.
(2) We numerically demonstrated the distinctive features of mirror diffusion and scattering diffusion.The diffusion of mirroring CRs is characterized by the oscillations in µ around µ = 0.As mirroring CRs repeatedly reverse their moving direction along the magnetic field, the mirror diffusion strongly enhances the confinement of CRs.By contrast, the scattering diffusion violates the first adiabatic invariant and is characterized by the stochastic change of µ within the range [µ c , 1].As scattering CRs do not reverse their moving direction, they are weakly confined and rapidly diffuse along the magnetic field.The transition between the two diffusion regimes occurs due to the pitch-angle scattering.
(3) We numerically measured the parallel mean free path of mirror diffusion, which can be much smaller than the turbulence injection scale and increases wtih increasing µ.
Our numerical measurements are in general consistent with the analytical predictions by LX21.As a new diffusion mechanism of CRs, the mirror diffusion has important implications.Around CR sources, e.g., supernova remnants (Xu (2021)), pulsar wind nebulae, the mirror diffusion of CRs with sufficiently large pitch angles can account for their slow diffusion, as indicated by gamma-ray observations (e.g., Abeysekara et al. (2017); Evoli et al. (2018)).In addition, the mirror diffusion of pickup ions is important for interpreting the Interstellar Boundary Explorer (IBEX) ribbon (Xu & Li (2023)).In the Galactic diffuse medium away from CR sources, CRs stochastically undergo slow mirror diffusion and fast scattering diffusion.This Lévy-flight-like propagation of CRs that we find is different from existing models of CR diffusion solely based on pitch-angle scattering.Its confrontation with CR observations will be further investigated.

Figure 2 .
Figure 2. Trajectory of a test particle in M2 with r L = 0.01Linj and µ0 = 0.05.The trajectory is color coded by magnetic field strength in (a) and time in (b).Trajectories in (a) and (b) are the same except that the trajectory in (a) is averaged over the gyroperiod Ω −1 to better illustrate the motion of the test particle.Arrows in (a) indicate the directions of mirror forces.

Figure 3 .
Figure 3. Gyro-averaged trajectory of a test particle in M2 with rL = 0.02Linj and µ0 = 0.05.The trajectory is color coded by time, with a transition from mirror diffusion to scattering diffusion.

Figure 4 .
Figure4.Time evolution of µ, particle position in z direction (the direction of the mean magnetic field for M2) z/Linj, and 2M B0/γmu 2 corresponding to the test particle trajectory in Fig.3.Blue and orange lines represent the numerical measurements and their gyro-averaged results, respectively.

Figure 5 .
Figure 5. Distributions of µ measured at different times for 5000 test particles with rL = 0.01Linj and µ0 = 0.05 (indicated by the vertical dotted line) in M2.The bin size of µ is 0.05.

Table 1 .
Parameters of MHD turbulence simulations and CRs are expected to undergo the mirror diffusion.At µ > µ c , scattering becomes dominant, and CRs are expected to undergo the scattering diffusion.