© The Authors 2025

1. Introduction

Cosmic rays (CRs), composed of high-energy charged particles, play an important role in astrophysics. Currently, the origin of Galactic CRs, particularly PeV CRs, is still an open issue (e.g., Amato & Casanova 2021). From an observational point of view, detecting the fluxes of secondary and primary CR particles can constrain the CR diffusion behavior (e.g., Amato 2014). Interestingly, Abeysekara et al. (2017) found that the CR diffusion coefficient in TeV halos surrounding Geminga and PSR B0656+14 is smaller by a factor of about 100 than the values constrained in the model of electron diffusion into the local interstellar medium (ISM). Studying the diffusion and acceleration of CRs in magnetized turbulent ISM is an effective way to understand the origin of the Galactic CRs. Meanwhile, it can help us understand some key astrophysical processes, such as the formation and evolution of galaxies and stars, the radiation of galaxy clusters, and diffuse γ-ray emissions (Rodgers-Lee et al. 2020; Semenov et al. 2021; Brunetti & Lazarian 2007; Krumholz et al. 2020; Yan et al. 2012).

Properly understanding the properties of the magnetohydrodynamic (MHD) turbulence is crucial for exploring the CR propagation processes. Traditionally, the slab/2D composite model of MHD turbulence was used to explain solar wind turbulence (e.g., Bieber et al. 1988; Matthaeus et al. 1990). Significantly, the new paradigm of MHD turbulence theory (Goldreich & Sridhar 1995, hereafter GS95) has advanced the study of the interaction of CRs with magnetized ISM (see Beresnyak & Lazarian 2019 for more details). Under the condition of the critical equilibrium, i.e., the time of the propagation of waves equal to the turbulence cascade time, the GS95 theory proposed that the motions of turbulence in the perpendicular and the parallel directions with respect to the magnetic field are different. The anisotropy of MHD turbulence has been verified by numerical simulations (Cho & Lazarian 2002, 2003, hereafter CL02, CL03). By decomposing MHD turbulence into three MHD modes: Alfvén, slow, and fast modes, CL02 found that the incompressible Alfvén mode and compressible slow mode are anisotropic, and the compressible fast mode presents isotropy.

Under the assumption of infinitesimal magnetic fluctuations, the quasi-linear theory (QLT) (Jokipii 1966) can describe the transport of particles. The QLT analytically predicted the relationship of λ_∥∝ρ^2−s (where s is the spectral index of the magnetic turbulence) between the parallel mean free path (MFP) and the rigidity ρ of CRs, which formulates a close link between CR transport behavior and magnetic turbulence properties (Hussein et al. 2015). However, when confronted with a realistic astrophysical environment with strong nonlinear turbulence, the QLT is limited in application. The most obvious drawback of the QLT is that it cannot deal with the 90° problem, resulting in theoretically infinite MFPs. To improve the theoretical description of the CR transport in realistic turbulent magnetic fields, various nonlinear theories (NLTs) and the extended QLT have been formulated (e.g., Voelk 1975; Matthaeus et al. 2003; Shalchi 2009; Yan & Lazarian 2008; Xu & Lazarian 2018). In the framework of the NLT, Qin (2002) numerically proposed the relation of λ_∥∝ρ^0.6, which is different from λ_∥∝ρ^1/3 of the QLT prediction for the Kolmogorov spectrum.

Given that CR diffusion depends on turbulent magnetic fields with anisotropy, it is reasonable to consider the perpendicular and parallel transport of CRs with regard to the (local) mean magnetic field. For the perpendicular diffusion measured with respect to the local mean magnetic field, it has been demonstrated that CRs undergo subdiffusion and superdiffusion in the dissipation range and inertial range, respectively (Hu et al. 2022; Zhang & Xu 2023). In the inertial range, the separation of particles increases with the time as t^3/2 due to the perpendicular superdiffusion of magnetic field lines (Lazarian & Yan 2014; Hu et al. 2022; Sampson et al. 2023). When considering the CR parallel diffusion, gyroresonance scattering, which describes the interaction of CRs with MHD turbulence (Yan & Lazarian 2008; Xu & Lazarian 2018), leads to scattering diffusion in the direction parallel to the local magnetic field. It has been claimed that the fast mode is a dominant scattering agent, and that slow and Alfvén modes with scale-dependent anisotropy can suppress the CR scattering (Yan & Lazarian 2002).

Except for the pitch-angle scattering by gyroresonant interactions, CRs are also reflected by magnetic mirrors (Fermi 1949; Cesarsky & Kulsrud 1973). Recently, Xu & Lazarian (2020) analytically predicted the mirroring and scattering effect of CRs in MHD turbulence. In addition, they found that the isotropic fast mode dominates both the mirroring and gyroresonant scattering of CRs in the case of compressible MHD turbulence, and the mirroring by slow mode dominates over the gyroresonant scattering by Alfvén and slow modes in the incompressible case. Given the intrinsic superdiffusion of the turbulent magnetic fields, which was predicted by Lazarian & Vishniac (1999), and numerically confirmed by Beresnyak (2013) with a high-resolution simulation (see also Lazarian et al. 2004), Lazarian & Xu (2021) proposed that CRs cannot be permanently trapped within one magnetic mirror, but move from one magnetic mirror to another along the magnetic field line. This is called mirror diffusion, which is slower than scattering diffusion. In addition, the most obvious advantage of mirror diffusion is that it can solve the 90° scattering problem.

Numerically, it has been demonstrated that CRs experience mirror diffusion at relatively large pitch angles and gyroresonant scattering at smaller pitch angles, the pitch-angle-dependent parallel MFP due to mirror diffusion alone being much smaller than the driving scale of MHD turbulence (Zhang & Xu 2023). At the same time, it has been found that lower-energy CRs preferentially undergo the mirror and wandering diffusion in the strong-field regions amplified by the nonlinear turbulent dynamo (Zhang & Xu 2024). Based on MHD simulations of star-forming regions, Barreto-Mota et al. (2024) also demonstrated that the parallel mirror diffusion dominates the CR diffusion behavior.

With the motivations of the slow diffusion found in observations and the theoretical predictions of turbulent mirror diffusion, we further numerically investigated the mirror and scattering diffusion behavior. One purpose of our work is to investigate the potential quantitative relationships between the CR MFPs and their energies. Another purpose is to study how individual MHD modes affect the mirror and scattering diffusion of CRs. The paper is organized as follows. Brief theoretical descriptions of scattering and mirroring are provided in Sect. 2. Section 3 describes our numerical methods, including the generation of 3D data of MHD turbulence, MHD mode decomposition, and the test-particle method. The numerical results are presented in Sect. 4, followed by a discussion and a summary in Sects. 5 and 6, respectively.

2. Theoretical description

Magnetohydrodynamic turbulence can be described by the ideal MHD equations (see also Sect. 3). In general, MHD turbulence can be decomposed into three modes: Alfvén, slow, and fast modes (e.g., Beresnyak 2019), numerically achieved by the Fourier transform (CL02) and the wavelet transform (Kowal & Lazarian 2010, hereafter KL10). It has been numerically demonstrated that Alfvén and slow modes have scale-dependent anisotropy, and the fast mode presents isotropy (CL02 and CL03). In addition, the Alfvén mode is incompressible, while the fast and slow modes are compressible.

The interactions of CR particles with turbulent magnetic fields affect the transport of CR particles in the directions parallel and perpendicular to the magnetic field. The resonant interaction of gyroresonance leads to scattering diffusion, while the nonresonant interaction with turbulent magnetic mirrors leads to mirror diffusion. Together they lead to parallel diffusion. According to the QLT theory, the linear resonance function of gyroresonance scattering is given by (Kulsrud & Pearce 1969; Voelk 1975)

$$R_n=\pi \delta (k_{\parallel }{v}_{\parallel }-\omega \pm n\mathrm{\Omega }),$$(1)

where ω and Ω are the wave frequency and particle gyrofrequency, respectively. This gyroresonance scattering requires that ω is Doppler-shifted to the gyrofrequency Ω of the particle or its cyclotron harmonics nΩ with a nonzero integer n (Yan & Lazarian 2008; Xu & Lazarian 2018). When the parallel velocity v_∥ of particles matches with the wave phase speed ω/k_∥, the resonance condition corresponds to a simplified resonance function

$$R_n=\pi \delta (k_{\parallel }{v}_{\parallel }-\omega ),$$(2)

resulting in efficient particle acceleration, i.e., the transit time damping (TTD) process (Yan & Lazarian 2002; Schlickeiser 2002). It is worth noting that TTD is one of the mechanisms for particle acceleration (Yan 2015).

In the case of nonresonant interaction, CRs can interact with the multi-scale mirrors naturally generated due to the turbulent compressions of magnetic fields. CR particles undergoing mirror diffusion satisfy the conditions of R_g<l_mir and μ<μ_c (Lazarian & Xu 2021), where R_g is the Larmor radius, μ the cosine of the pitch angle, and l_mir the scale of the mirror. The pitch angle is the angle between the velocity of the particle and the magnetic field. We note that the critical value μ_c can be obtained by μ_c=min[μ_mir,μ_eq], where ${\mu }_{\mathrm{mir}}\approx \sqrt{\frac{\delta B}{B_0+\delta B}}\approx \sqrt{\frac{\delta B}{B_0}}$ corresponds to the loss cone angle (B₀ and δB denote the mean and fluctuating magnetic fields, respectively), and ${\mu }_{\mathrm{eq}}\simeq [\frac{14}{\pi }(\frac{\delta B_{\mathrm{f}}}{B_0}{)}^2(\frac{R_{\mathrm{g}}}{L_{\mathrm{inj}}}{)}^{\frac{1}{2}}{]}^{\frac{2}{11}}$ corresponds to μ value when the mirroring rate equals to scattering rate of particles. In the case of compressible MHD turbulence, due to the fast mode dominating the mirror diffusion (i.e., the mirroring rate of CRs by fast mode greater than that by slow mode), μ_mir can be replaced by (Lazarian & Xu 2021)

$${\mu }_{\mathrm{mir}}=\sqrt{\frac{\delta B_{\mathrm{f}}}{B_0+\delta B_{\mathrm{f}}}}·$$(3)

The mirroring particles preserve the first adiabatic invariant, that is, the conservation of the magnetic moment $M=\frac{\gamma {\mathit{mu}}_{\perp }^2}{2B}=\mathrm{const}$, where γ, m, and u_⊥ denote the Lorentz factor, mass, and perpendicular velocity of CR particles, respectively.

The perpendicular transport of CRs is influenced by the superdiffusion of the magnetic field. Since the perpendicular diffusion of CRs depends on the transport scale, one can divide the perpendicular transport process into two parts: the large-scale transport (the scale larger than the injection scale L_inj), and the small-scale transport (the scale smaller than L_inj). According to the M_A values (see Sect. 3 for its definition), the turbulence can be divided into super-Alfvénic (M_A>1) and sub-Alfvénic (M_A<1) turbulence regimes. For the former the transition from hydrodynamic turbulence to MHD turbulence happens at $L_{\mathrm{A}}=L_{\mathrm{inj}}{M}_{\mathrm{A}}^{-3}$, while for the latter the transition from weak turbulence to strong turbulence happens at $L_{\mathrm{tr}}=L_{\mathrm{inj}}{M}_{\mathrm{A}}^2$ (see Beresnyak & Lazarian 2019 for more details).

The propagation of CRs is related to the magnetization degree of turbulence. In the case of super-Alfénic turbulence, on a scale larger than L_A, with turbulent energy dominating the magnetic energy, the diffusion is expected to be isotropic, which means that the parallel diffusion coefficient D_∥ and the perpendicular diffusion coefficient D_⊥ are similar (Yan & Lazarian 2008; Lazarian & Yan 2014):

$$D_{\perp }\approx D_{\parallel }.$$(4)

However, when the scale is less than L_A, the diffusion is expected to be anisotropic, and the parallel and perpendicular diffusion coefficients follow (Maiti et al. 2022)

$$D_{\perp }\approx D_{\parallel }{M}_{\mathrm{A}}^3.$$(5)

For the sub-Alfvénic turbulence, the parallel and perpendicular diffusion coefficients maintain the relation (Lazarian & Yan 2014)

$$D_{\perp }\approx D_{\parallel }{M}_{\mathrm{A}}^4,$$(6)

which has been numerically confirmed by Xu & Yan (2013) and Maiti et al. (2022).

3. Simulation methods

To model turbulent magnetic fields, we performed 3D numerical simulations to solve the following ideal MHD equations:

$$\frac{\partial \rho }{\partial t}+\nabla ·(\rho v)=0,$$(7)

$$\rho \left[\frac{\partial v}{\partial t}+(v·\nabla )v\right]+\nabla p-\frac{1}{4\pi }J\times B=f,$$(8)

$$\frac{\partial B}{\partial t}-\nabla \times (v\times B)=0,$$(9)

$$\nabla ·B=0.$$(10)

Here f is a random force to drive the turbulence, J=∇×B the current density, and v the fluid velocity. In Eq. (8), $p=c_{\mathrm{s}}^2\rho$ represents gas pressure, where c_s and ρ are the sonic speed and density, respectively.

We use a second-order-accurate hybrid simulation code developed by Cho & Lazarian (2002), which is essentially non-oscillatory. Considering a periodic boundary condition and setting a non-zero mean magnetic field in the x-axis direction, we drive 3D simulations by a pure solenoidal forcing in the Fourier space. With the numerical resolution of L³ = 512³, we set the injection wave number k_in≃2.5, corresponding to the injection scale of 0.4L. When the simulation reaches a statistically steady state (with a steady turbulent energy spectrum), we take 3D data cubes with density, velocity, and magnetic field. Defining the Alfvénic Mach number $M_{\mathrm{A}}=\langle \frac{|v|}{V_{\mathrm{A}}}\rangle$ and the sonic Mach number $M_{\mathrm{s}}=\langle \frac{|v|}{c_{\mathrm{s}}}\rangle$, we can obtain the values of the characteristic parameters (listed in Table 1) of each simulation. Here, the Alfvén speed V_A is related to the total magnetic field strength and plasma density by $V_{\mathrm{A}}=\frac{|B|}{\sqrt{\rho }}$. In addition, the plasma parameter $\beta =p_{\mathrm{gas}}/p_{\mathrm{mag}}=2M_{\mathrm{A}}^2/M_{\mathrm{s}}^2$ is defined as the ratio of gas pressure to magnetic pressure.

We decompose the compressible MHD turbulence into Alfvén, slow, and fast modes by the wavelet transforms (KL10). The displacement vectors of three modes are defined as (CL02; KL10)

$$\widehat{{\xi }_{\mathrm{f}}}\propto (1+\alpha +\sqrt{D})k_{\perp }\widehat{k_{\perp }}+(-1+\alpha +\sqrt{D})k_{\parallel }\widehat{k_{\parallel }},$$(11)

$$\widehat{{\xi }_{\mathrm{s}}}\propto (1+\alpha -\sqrt{D})k_{\perp }\widehat{k_{\perp }}+(-1+\alpha -\sqrt{D})k_{\parallel }\widehat{k_{\parallel }},$$(12)

$$\widehat{{\xi }_{\mathrm{A}}}\propto \widehat{k_{\perp }}\times \widehat{k_{\parallel }},$$(13)

where $D=(1+\alpha {)}^2-4\alpha {cos}^2\theta ,\alpha =c_{\mathrm{s}}^2/V_{\mathrm{A}}^2$ and θ is the angle between the wave vector k and B₀. We first use the discrete wavelet transform to obtain the wavelet coefficients corresponding to the magnetic field. Next, we perform the Fourier decomposition for three modes following CL02 by projecting the magnetic field into the displacement vectors of three modes (see Eqs. (11)–(13)) in the Fourier space. Finally, we carry out the inverse Fourier transform for three modes to get the information on the magnetic field in real space.

By injecting test particles into 3D data cubes, we can simulate the transport of CR particles. Numerically, we use the Bulirsch–Stoer method (Press et al. 1986) with an adaptive time step to trace the motions of the charged particles. Specifically, the motions of the test particles are governed by the Lorentz equation as

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{mu}{\sqrt{(1-u^2/c^2)}}\right)=q(u\times B),$$(14)

where q and $\gamma =\frac{1}{\sqrt{(1-u^2/c^2)}}$ represent the charge and Lorentz factor of the particle, respectively. The particle position r is related to its velocity by

$$\frac{\mathrm{d}r}{\mathrm{d}t}=u,$$(15)

where the time t in code units is normalized as τ=t·Ω and Ω=u/R_g. In our simulation, we use the ratio of the Larmor radius R_g to L_inj to characterize test particle energy, where R_g is defined as

$$R_{\mathrm{g}}=\gamma \frac{m{c}^2}{q{B}_0}·$$(16)

To numerically describe the diffusion processes of CR particles, we define the parallel diffusion coefficient D_∥ and the perpendicular coefficient D_⊥ with respect to the mean magnetic field direction as (Giacalone & Jokipii 1999)

$$D_{\parallel }=\frac{\langle \mathrm{\Delta }{\ell }_{\parallel }^2\rangle }{2\mathrm{\Delta }t}$$(17)

and

$$D_{\perp }=\frac{\langle \mathrm{\Delta }{\ell }_{\perp }^2\rangle }{2\mathrm{\Delta }t},$$(18)

respectively. In Eqs. (17) and (18), the particle's parallel and perpendicular displacements squared with respect to the mean magnetic field are defined as

$$\mathrm{\Delta }{\ell }_{\parallel }^2=(x-x_0{)}^2$$(19)

and

$$\mathrm{\Delta }{\ell }_{\perp }^2=(y-y_0{)}^2+(z-z_0{)}^2,$$(20)

respectively. The symbol 〈…〉 denotes the average over all particles. Furthermore, we can define the parallel and perpendicular MFPs of particles as

$${\lambda }_{\parallel }=\frac{3D_{\parallel }}{u}\hspace{0.5em}\hspace{0.5em}\mathrm{and}\hspace{0.5em}\hspace{0.5em}{\lambda }_{\perp }=\frac{3D_{\perp }}{u},$$(21)

respectively.

4. Numerical results

4.1. Spatial diffusion of single particle

Before exploring the statistical properties of particle diffusion, here we focus on the behaviors of the individual initially mirroring and scattering particles. Figure 1 depicts the trajectories of the randomly selected initially mirroring (left panel) and initially scattering (right panel) particles. We would like to note that the mirroring and scattering particles are prepared initially, but at later times, they will lose their identification as μ changes through scattering. Due to the setting of the mean magnetic field along the x-axis direction, we see that the curved magnetic field lines (thin solid lines) distribute along the x-axis. We set the initial pitch-angle cosines μ₀ = 0.15 and μ₀ = 0.8 corresponding to the initially mirroring and scattering particles, respectively. Their trajectories are color-coded by time τ. Compared to these two particles in Fig. 1, we can see that the parallel spatial extent of the initially mirroring particle (left panel) is smaller than that of the initially scattering particle (right panel). Physically, the mirroring effect can more effectively confine particles than scattering can. Due to the perpendicular superdiffusion of magnetic field lines, mirroring particles can interact with different magnetic mirrors by bouncing back and forth among several different magnetic mirrors. At a later time, the initially mirroring particle changes its pitch angle due to the scattering. When the pitch angle decreases to a sufficiently small value, i.e., close to the condition of μ>μ_c, the mirroring particle escapes from turbulent magnetic mirrors and undergoes scattering diffusion. When the pitch angle becomes sufficiently large with μ<μ_c, it will again undergo mirror diffusion.

Fig. 1. Trajectories of the initially mirroring particle (left panel; initial pitch-angle cosine μ0 = 0.15) and the initially scattering particle (right panel; μ0 = 0.8). The magnetic field lines (black) extend to 1024 pixels in the x-axis direction. The trajectories of these two particles are color-coded by the time τ (normalized by the gyrofrequency Ω), as shown in the color bar. We note that for the sake of comparison, the right panel only shows a part of the trajectory in the range of τ<160. The Larmor radius of these two particles is Rg = 0.03Linj. The simulations are from R2 listed in Table 1.

To observe the detailed features of these two particles, we plot μ, magnetic moment M, and the position in the x-axis direction L_x as a function of time in Fig. 2. Panels a and b correspond respectively to the initially mirroring and scattering particles in Fig. 1. As seen in panel a, the initially mirroring particle undergoes mirror diffusion before τ≃220 with an insignificant change in L_x and the pitch angle crossing 90° over and over again. When μ approaches the critical condition μ_c, the particle undergoes a weak scattering process, making it escape from the current mirror and enter another mirror. Again, the particle experiences a similar mirror diffusion in the range from τ≃300 to 700. We note that during the mirror diffusion process, the (normalized) magnetic moment maintains M≃1.0. After τ≃750, the particle experiences a significant scattering that allows the particle to leave the mirrors; it then undergoes scattering diffusion. If one traces the particle long enough, it can undergo mirror diffusion again when the conditions of mirror diffusion: μ<μ_c and M=const are satisfied. It is found that the combination of the mirror and scattering diffusion results in a Lévy-flight-like propagation for CRs (see also Barreto-Mota et al. 2024 for more details).

Fig. 2. Cosine of pitch angle μ, normalized magnetic moment 2MB0/γmu2, and spatial displacement Lx/Linj in the x-axis direction as a function of time τ. The results of the initially mirroring and scattering particles shown in panels a and b correspond respectively to the left and right particles in Fig. 1. The horizontal dashed lines represent ±μc in Table 1. The blue and red curves represent respectively the numerical result and the averaged result within one gyration period.

As seen in panel b, the initially scattering particle rapidly diffuses a long distance along the magnetic field before τ≃300, while μ and M almost maintain their initial values, i.e., μ≃0.8 and M≃0.3, respectively. This indicates inefficient scattering, and thus the particle simply travels along field lines. It is worth noting that for the initially mirroring and scattering particles, they both remain constant in the magnetic moment at the early stage of time, but M≃0.3 for the initially scattering particle is different from M≃1.0 for the initially mirroring particle. Therefore, one cannot use the constant magnetic moment alone to identify a mirroring particle¹. In the range from τ≃400 to 700, the particle satisfies the conditions of mirror diffusion mentioned above, with the insignificant change in L_x. After τ≃700, the particle briefly undergoes scattering diffusion and then goes back to the mirror diffusion region.

4.2. Parallel and perpendicular displacements of particles

This section explores statistically the diffusion properties of a large sample of initially mirroring and scattering particles. Figure 3 shows the results of the parallel displacement $\langle \mathrm{\Delta }{\ell }_{\parallel }^2\rangle$ (panel a) and the perpendicular displacement $\langle \mathrm{\Delta }{\ell }_{\perp }^2\rangle$ (panel b) of the initially mirroring and scattering particles as a function of time, and the ratio of $\langle \mathrm{\Delta }{\ell }_{\parallel }^2\rangle$ to $\langle \mathrm{\Delta }{\ell }_{\perp }^2\rangle$ (panels c and d) for them, using the models on R1, R2, and R3 listed in Table 1. We set μ₀∈[0.1,0.2] and μ₀∈[0.8,0.9] respectively for the initially mirroring and scattering particles with the same Larmor radius R_g = 0.03L_inj. In both cases, we injected 2000 test particles to trace their trajectories. After the particles underwent thousands of gyroperiods, they reached a normal diffusion regime². We terminated our test particle simulation.

Fig. 3. Upper panels: Mean square displacement traveled by overall initially mirroring and scattering particles vs. the time (normalized by the gyrofrequency Ω) in the parallel (panel a) and perpendicular (panel b) directions with respect to the mean magnetic fields. Lower panels: Ratio of the parallel (panel c) to perpendicular (panel d) mean square displacement between the initially scattering and mirroring particles. The initial pitch angle values are set to μ0∈[0.1,0.2] and μ0∈[0.8,0.9] for initially mirroring and scattering particles, respectively. The Larmor radius of particles is Rg = 0.03Linj.

In the case of parallel diffusion (see panels a and c), the particles’ diffusion can be divided into two stages in terms of the dependence of mean square displacement on time. The first stage is superdiffusion with a power-law index larger than 1. The second stage is normal diffusion, where the mean square displacement presents a linear dependence on time of $\langle \mathrm{\Delta }{\ell }_{\parallel }^2\rangle \propto \tau$. In the first stage, we can see $\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{m}}^2\rangle <\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{s}}^2\rangle$ for three cases from panel c, which means the diffusive displacement of the initially mirroring particles is smaller than that of the initially scattering particles during the same time. This is because the initially mirroring particles are more effectively confined compared to the initially scattering particles. It means that the presence of magnetic mirrors makes the mirror diffusion confine CRs more strongly than scattering diffusion. In the second stage, we can see $\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{m}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{s}}^2\rangle$. It indicates that their displacements converge, irrespective of the initial pitch angles, when the parallel diffusion is undergoing normal diffusion. Compared with different turbulence regimes, the main factor affecting parallel diffusion is M_A rather than M_s. As seen in the first stage, the initially mirroring and scattering particles have a larger number of gyrations and faster parallel diffusion for sub-Alfvénic turbulence (see R2 and R3 listed in Table 1) than that for super-Alfvénic turbulence (see R1). As for the case of R1, due to the enhanced fluctuation of the magnetic field, particles are subjected to more efficient scattering of the magnetic field, the initially mirroring and scattering particles earlier to reach normal diffusion.

In the case of perpendicular diffusion (see panels b and d), except for the presence of a significant oscillation at τ≲5, which may be from numerical dissipation, we can also see $\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{m}}^2\rangle <\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{s}}^2\rangle$ for superdiffusion and $\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{m}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{s}}^2\rangle$ for normal diffusion in the different turbulence regimes. Compared with R1, R2, and R3 within the same time, the initially mirroring and scattering particles all diffuse larger distances for R1 than those for R2 and R3, for which the distance is comparable. In addition, we note that R1 shows the isotropic diffusion of $\langle \mathrm{\Delta }{\ell }_{\perp }^2\rangle \approx \langle \mathrm{\Delta }{\ell }_{\parallel }^2\rangle$ during the whole process, whereas R2 and R3 indicate anisotropic diffusion of $\langle \mathrm{\Delta }{\ell }_{\perp }^2\rangle <\langle \mathrm{\Delta }{\ell }_{\parallel }^2\rangle$. The reason is that the perpendicular diffusion of particles is subject to that of magnetic fields. With increasing the fluctuation of the magnetic field, as characterized via δB_rms/〈B〉 listed in Table 1, the diffusion of CRs would exhibit isotropic property, which is consistent with the expected isotropic diffusion of CRs for super-Alfvénic turbulence.

To further understand the effect of scattering on the distribution of particle pitch angles, we plot the distribution of μ corresponding to the initially mirroring particles and scattering particles at different times during the diffusion process. Using model R2 listed in Table 1, we plot the probability density function of μ in Fig. 4 corresponding to Fig. 3. As shown, at the first gyroperiod, i.e., τ = 1, the peak of μ happens at μ≃0.0 for initially mirroring particles (see panel a) and μ≃0.8 for initially scattering particles (see panel b). These two peaks approximately approach the initial μ distribution that we set. When increasing the time to τ = 10⁴, we see that μ is evenly distributed between −1 and 1 for the initially mirroring and scattering particles, that is, the distribution of μ is gradually randomized. This demonstrates that, due to scattering, irrespective of the initial pitch angle, particles stochastically undergo mirroring and scattering in turn, resulting in similar diffusion behavior over a sufficiently long time.

Fig. 4. Probability density function of μ at different times τ, arising from R2 in Fig. 3. The vertical dashed lines correspond to μ = 0 in panel a and μ = 0.8 in panel b.

4.3. Mean free path of particles

In this section we investigate the dependence of the propagation of initially mirroring and scattering particles on their energies in different turbulence regimes. Figure 5 plots the parallel and perpendicular MFPs, λ_∥ and λ_⊥ (see Eq. (21)), as a function of the Larmor radius, R_g, where the error bars plotted represent the standard deviation of the MFPs. With the same initial setup of μ as in Fig. 3, we ran a series of simulations by changing R_g values from the dissipation scale (several grids) to the injection scale L_inj. In the early stage (around several hundred gyroperiods), we see superdiffusion for all the simulations. During the middle and later stages, where the initially mirroring and scattering particles reach the normal diffusion, i.e., 〈Δℓ²〉∝τ, we measured the parallel and perpendicular MFPs via the parallel and perpendicular diffusion coefficient (see also Sect. 3).

Fig. 5. MFPs of mirroring and scattering particles measured at different CR energies. The initial pitch angle values are set to μ0∈[0.1,0.2] and μ0∈[0.8,0.9] for mirroring and scattering particles, respectively. The results in panels a, b, and c are based on the R1, R2, and R3 listed in Table 1, respectively. The different symbols denote the numerical results, and the solid lines are a linear fitting for them. The vertical dash-dotted lines indicate the transition scales. The blue area plotted in each panel represents the dissipation scales estimated by the power spectra of the magnetic field.

In the case of the super-Alfvénic turbulence, the numerical result is shown in Fig. 5a, from which we see that the relation between λ and R_g presents two significant features separated by the scale L_A describing the transition from the hydrodynamic turbulence to strong MHD turbulence (see Sect. 2). In the range of L_dis<R_g<L_A, the linear fitting for parallel and perpendicular MFPs versus R_g shows the relation of ${\lambda }_{\parallel }\simeq {\lambda }_{\perp }=5.12R_{\mathrm{g}}^{1.04}\propto R_{\mathrm{g}}$. Similarly, we can also see that in the dissipation region of R_g<L_dis, the parallel and perpendicular MFPs are approximately proportional to R_g. Given that the random magnetic field dominates the regular magnetic field, δB_rms/〈B〉≃1.08, in this simulation, a strong random fluctuation of the magnetic field lines leads to the absence of anisotropic diffusion.

In the high-energy range of R_g>L_A, we observe a significant difference between the parallel and perpendicular MFPs. The former satisfies a steeper power-law relationship of ${\lambda }_{\parallel }=16.59R_{\mathrm{g}}^{1.85}\propto R_{\mathrm{g}}^2$, while the latter reaches a plateau-like distribution with increasing energy. Here, the turbulence cascade essentially retains its fluid-like properties. The plateau-like distribution of λ_⊥ and the steeper power-law relation of λ_∥ versus R_g may be due to the weak magnetic field having a marginal influence on turbulence. In addition, the plateau-like distribution of λ_⊥ versus R_g can also be associated with the magnetic field wandering. Compared with the initially scattering and mirroring particles with the same energy, we note that the parallel or perpendicular MFPs maintain the same value. It indicates that the initial μ of particles does not change the MFPs of particles.

As for the sub-Alfvénic turbulence, our results are plotted in Figs. 5b and c. As is shown, the perpendicular MFPs present a power-law distribution of ${\lambda }_{\perp }\propto R_{\mathrm{g}}^{0.7}\propto R_{\mathrm{g}}^{2/3}$ in the low-energy range, following a plateau beyond the transition scale L_tr. This power-law relation is slightly shallower than λ_⊥∝R_g shown in Fig. 5a. Generally, the power-law relation of the perpendicular diffusion maintains a similarity between sub-Alfvénic and super-Alfvénic turbulence. The significant difference is from the parallel diffusion of particles. First, the relation between λ_∥ and R_g presents an inverse distribution in the range of the dissipation scale, i.e., the negative index of λ_∥ versus R_g. This may be associated with the propagation of CRs in damped turbulence (e.g., Yan & Lazarian 2008). Second, the linear fitting provides a power-law relation of $\lambda \propto R_{\mathrm{g}}^{0.35}\propto R_{\mathrm{g}}^{1/3}$ in the range of strong turbulence of L_dis<R_g<L_tr. Third, the diffusion of the initially mirroring and scattering particles exhibits a distinctly anisotropic characteristic, with the parallel MFPs greater by two to three orders of magnitude than the perpendicular ones at the same CR energy.

Similarly to Fig. 5a in the range of R_g>L_tr, Figs. 5b and c present a steep power-law of ${\lambda }_{\parallel }\propto R_{\mathrm{g}}^{ϵ}$, with ϵ≃1.89 (panel b) and 1.63 (panel c). The plateau distribution of λ_⊥ demonstrates that in the weak turbulence regime, the parallel diffusion of initially scattering and mirroring particles behaves similarly due to the wandering of the field lines. As a result, we find that the perpendicular and parallel diffusion of initially mirroring and scattering particles strongly depend on the magnetized parameter M_A of MHD turbulence.

4.4. Influence of MHD modes on spatial diffusion

To explore the influence of different MHD modes on diffusion processes of the initially mirroring and scattering particles, we first decomposed 3D MHD turbulence data (R1 as an example) using the wavelet decomposition method (see Sect. 3) and then injected test particles into the post-decomposed data of individual modes in the same procedure mentioned above. As shown in Fig. 6, the Alfvén and slow modes present the anisotropy scaling of ${\mathcal{L}}_{\parallel }\propto {\mathcal{L}}_{\perp }^{2/3}$, with the power spectrum of E∝k^−5/3. The fast mode has an isotropy scaling of ℒ_∥∝ℒ_⊥, with the index of power spectrum close to −2 in the inertial range due to the shockwave. In addition, the magnetic energies of three modes satisfy the following relationship E_B,S>E_B,A>E_B,F.

Fig. 6. Power spectra (panel a) and anisotropy scalings (panel b) of magnetic fields corresponding to the Alfvén, fast, and slow modes decomposed from R1.

Figures 7a and b show the evolution of the parallel and perpendicular mean square displacements of the initially mirroring and scattering particles. As shown, the parallel and perpendicular diffusion of particles in three modes experience the transition from superdiffusion to normal diffusion, similar to those in the pre-decomposed case. The magnetosonic modes dominate the parallel diffusion processes for the initially mirroring and scattering particles (see panel a), and the Alfvén mode dominates the perpendicular diffusion processes (see panel b). Compared with initially mirroring and scattering particles, we note that significant differences occur in the early stages of time (τ<10²).

Fig. 7. Upper panels: Mean square displacement of the parallel (left) and perpendicular (right) diffusion of the initially mirroring and scattering particles vs. the time (normalized by Ω) for the Alfvén (A), fast (F), and slow (S) modes. Lower panels: Ratio of the parallel (left) to perpendicular (right) mean square displacement between their three plasma modes. The subscripts i and j represent any two of Alfvén, fast, and slow modes. The simulations with the Larmor radius of Rg = 0.03Linj are based on the decomposed R1 data. The μ values for initially mirroring and scattering particles are set to μ0∈[0.1,0.2] and μ0∈[0.8,0.9], respectively.

Moreover, in Figs. 7c and d, we plot the ratio of the parallel and perpendicular mean square displacements between different modes, respectively. In the case of parallel diffusion (panel c), we can see $\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{F}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{S}}^2\rangle >\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{A}}^2\rangle$ for the initially scattering particles during the whole period. However, for the initially mirroring particles, we see $\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{F}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{S}}^2\rangle <\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{A}}^2\rangle$ at τ≲6, and $\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{F}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{S}}^2\rangle >\langle \mathrm{\Delta }{\ell }_{\parallel ,\mathrm{A}}^2\rangle$ at τ≳6. In the case of perpendicular diffusion (panel d), we can see that the diffusion of initially mirroring particles satisfies $\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{A}}^2\rangle <\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{S}}^2\rangle <\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{F}}^2\rangle$ at τ≲10 and $\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{A}}^2\rangle >\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{S}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{F}}^2\rangle$ at τ≳10. As for the initially scattering particles, we can see $\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{A}}^2\rangle >\langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{S}}^2\rangle \simeq \langle \mathrm{\Delta }{\ell }_{\perp ,\mathrm{F}}^2\rangle$ during the whole periods.

In short, the magnetosonic and the Alfvén modes dominate the parallel and perpendicular diffusion respectively of the mirroring and scattering particles for a long evolution period. In the magnetosonic modes the parallel diffusion of particles is faster, whereas the perpendicular diffusion of particles is faster in the Alfvén mode. Parallel diffusion is associated with the compressibility of the slow and fast modes, while perpendicular diffusion is associated with the random walk of the magnetic field dominated by the Alfvén mode.

5. Discussion

In this work, we numerically explored the mirror and scattering diffusion of CRs in sub-Alfvénic and super-Alfvénic turbulence regimes first. After particles reached normal diffusion within several thousands of gyroperiods in MHD turbulence, we measured the parallel and perpendicular MFPs of the CRs and quantified the relationships between MFPs and the Larmor radii. We found that the particles with larger initial pitch angles (μ₀∈[0.1,0.2]) diffuse more slowly than the particles with smaller initial pitch angles (μ₀∈[0.8,0.9]), which is consistent with Barreto-Mota et al. (2024). The MFP of CRs is independent of the initial pitch angles. After testing the influence of local and global reference frames on the measurement results, we found that, for the power-law stages, the change in the frame of reference does not affect the power-law relation we obtained.

As predicted theoretically in Lazarian & Yan (2014), the CR perpendicular superdiffusion with the relation of 〈Δℓ_⊥〉∝t^3/2 originates from the superdiffusion properties of turbulent magnetic field lines in the local frame of reference. Numerically, the perpendicular superdiffusion was verified by Xu & Yan (2013), Hu et al. (2022), and Maiti et al. (2022). We note that Maiti et al. (2022) explored CR superdiffusion behavior in the global and local references and found that the index of CR superdiffusion depends slightly on the choice of the frame of reference. Given that one cannot observe CR diffusion behavior in the local frame of reference, all the results presented in this work are based only on the global frame of reference. In general, the efficiency of resonant scattering decreases for the same M_A with an increase in the inertial range. The results presented in this paper are limited to the numerical resolution of 512³ (see Zhang & Xu 2023 for a comparison between different resolutions).

Our current studies do not consider the CR acceleration and its radiative losses when the mirroring and scattering effects happen. Cohet & Marcowith (2016) has studied the propagation of CRs in MHD turbulence considering different forcing effects. They also measured the dependence of the MFP of CRs on their energy. For the sub-Alfvénic regime, they found the relations of ${\lambda }_{\parallel }\propto R_{\mathrm{g}}^{1/3}$ for M_A>0.5 and ${\lambda }_{\perp }\propto R_{\mathrm{g}}^{0.67}$ for M_A = 0.67, which are consistent with the results presented in our paper (see the results for R2 and R3). Similarly, they also see inverted spectra with solenoidal forcing when M_A<1. In our work, we focus on the influence of different turbulence regimes on the propagation of CRs and measure MFPs of CRs within a wide energy range, i.e., from several grids to L_inj. Considering the effect of acceleration on particles propagation, Beresnyak et al. (2011) found that the parallel diffusion coefficient is proportional to R_g for R_g<0.1L_inj and $R_{\mathrm{g}}^2$ for L_inj>R_g>0.1L_inj in trans-Alfvénic incompressible MHD turbulence. In addition, they also found that the perpendicular diffusion coefficient is independent of energy.

Our current work also explored the interaction of CRs with three MHD modes during the mirror and scattering diffusion. We found that the magnetosonic modes dominate parallel diffusion and the Alfvén mode dominates perpendicular diffusion. In our previous works, we studied the CR acceleration, diffusion, and scattering of CRs due to the interactions of individual MHD modes. Zhang & Xiang (2021) found that the fast mode dominates the acceleration of particles in the case of super-Alfvénic and supersonic turbulence and the slow mode dominates the acceleration for sub-Alfvénic turbulence. As for the diffusion of the accelerated CRs, the slow mode dominates the diffusion of particles in the strong turbulence regime, whereas three modes have a comparable role in the weak turbulence regime (Gao & Zhang 2024).

In addition to the mirror diffusion explored in this paper, the understanding of the slow diffusion phenomena around the source region is attempted in other ways including the anisotropy diffusion of CRs due to the anisotropy of MHD turbulence (e.g., Liu et al. 2019; Gao & Zhang 2025); the self-generated turbulence caused by CR flow restricting the diffusion of CRs (e.g., Evoli et al. 2018); the constraint of diffusion from the strong turbulence generated by the shock wave of the parent supernova remnant (Fang et al. 2019); and the two-zone mode assuming different diffusion coefficient near and far from the source region (Fang et al. 2018). In this paper, we confirmed that CRs undergo not only nonresonant interaction of magnetic mirrors, but also gyroresonance interactions with the magnetic field during diffusion, exhibiting a synergistic effect of mirror and scattering diffusion (Barreto-Mota et al. 2024), which provides an alternative way to understand the diffusion process of CRs in different environments such as supernova remnants (see also Xu 2021 for application of mirror diffusion) and massive star-forming regions.

6. Summary

With the modern understanding of MHD turbulence, we performed test-particle simulations to study the mirror diffusion and scattering diffusion of CRs. The main findings of the paper are summarized as follows:

1. Our simulations demonstrate that CRs with large pitch angles undergo mirror diffusion due to mirror reflection. In contrast, CRs with small pitch angles experience scattering diffusion due to gyroresonance with the magnetic field. The transition between the mirror and scattering diffusion appears when meeting or violating the following condition: the magnetic mirror size is larger than the Larmor radius of the particle with μ<μ_c and M=const. A comparison of mirror and scattering diffusions shows that the former has a stronger confining effect on CRs.

2. Regardless of the setting of the initial pitch angles, we find that CRs experience a transition from a superdiffusion to a normal diffusion after a sufficiently long evolution. During superdiffusion, the mean square displacement corresponding to the initially mirroring particles is significantly smaller than that corresponding to the initially scattering particles.

3. The normal diffusion stage is involved in a mixture of scattering and mirror diffusion, where the interaction of CRs with MHD turbulence results in a significant anisotropy of CR diffusion in the sub-Alfvénic regime.

4. The CR diffusion strongly depends on the properties of MHD turbulence. We find the presence of the power-law relation between the mean free path and the Larmor radius R_g. When the R_g of CRs is in the range of strong turbulence, we have ${\lambda }_{\perp }\propto R_{\mathrm{g}}^{2/3}$ and ${\lambda }_{\parallel }\propto R_{\mathrm{g}}^{1/3}$ for the sub-Alfvénic regime as well as λ_⊥≃λ_∥∝R_g for the super-Alfvénic regime. When the R_g of CRs is in the range of weak turbulence (M_A<1) or hydrodynamic turbulence (M_A>1), we find the relation of ${\lambda }_{\parallel }\propto R_{\mathrm{g}}^2$, and the plateau-like distribution of λ_⊥.

5. Compressibility of the magnetosonic modes dominates the parallel diffusion of CR particles, while the random walk of the magnetic field associated with the Alfvén mode dominates the perpendicular diffusion of CR particles.

Acknowledgments

We thank the anonymous referee for valuable comments that significantly improved the quality of the paper. We thank Jungyeon Cho for the helpful discussions on the referee's comments. Y.W.X. thanks Chao Zhang for the helpful discussions on numerical methods of test-particle simulations. J.F.Z. is grateful for the support from the National Natural Science Foundation of China (No. 12473046) and the Hunan Natural Science Foundation for Distinguished Young Scholars (No. 2023JJ10039).

References

All Tables

All Figures

Fig. 1. Trajectories of the initially mirroring particle (left panel; initial pitch-angle cosine μ0 = 0.15) and the initially scattering particle (right panel; μ0 = 0.8). The magnetic field lines (black) extend to 1024 pixels in the x-axis direction. The trajectories of these two particles are color-coded by the time τ (normalized by the gyrofrequency Ω), as shown in the color bar. We note that for the sake of comparison, the right panel only shows a part of the trajectory in the range of τ<160. The Larmor radius of these two particles is Rg = 0.03Linj. The simulations are from R2 listed in Table 1.

In the text

Fig. 2. Cosine of pitch angle μ, normalized magnetic moment 2MB0/γmu2, and spatial displacement Lx/Linj in the x-axis direction as a function of time τ. The results of the initially mirroring and scattering particles shown in panels a and b correspond respectively to the left and right particles in Fig. 1. The horizontal dashed lines represent ±μc in Table 1. The blue and red curves represent respectively the numerical result and the averaged result within one gyration period.

In the text

Fig. 3. Upper panels: Mean square displacement traveled by overall initially mirroring and scattering particles vs. the time (normalized by the gyrofrequency Ω) in the parallel (panel a) and perpendicular (panel b) directions with respect to the mean magnetic fields. Lower panels: Ratio of the parallel (panel c) to perpendicular (panel d) mean square displacement between the initially scattering and mirroring particles. The initial pitch angle values are set to μ0∈[0.1,0.2] and μ0∈[0.8,0.9] for initially mirroring and scattering particles, respectively. The Larmor radius of particles is Rg = 0.03Linj.

In the text

	Fig. 4. Probability density function of μ at different times τ, arising from R2 in Fig. 3. The vertical dashed lines correspond to μ = 0 in panel a and μ = 0.8 in panel b.
In the text

Fig. 5. MFPs of mirroring and scattering particles measured at different CR energies. The initial pitch angle values are set to μ0∈[0.1,0.2] and μ0∈[0.8,0.9] for mirroring and scattering particles, respectively. The results in panels a, b, and c are based on the R1, R2, and R3 listed in Table 1, respectively. The different symbols denote the numerical results, and the solid lines are a linear fitting for them. The vertical dash-dotted lines indicate the transition scales. The blue area plotted in each panel represents the dissipation scales estimated by the power spectra of the magnetic field.

In the text

	Fig. 6. Power spectra (panel a) and anisotropy scalings (panel b) of magnetic fields corresponding to the Alfvén, fast, and slow modes decomposed from R1.
In the text

Fig. 7. Upper panels: Mean square displacement of the parallel (left) and perpendicular (right) diffusion of the initially mirroring and scattering particles vs. the time (normalized by Ω) for the Alfvén (A), fast (F), and slow (S) modes. Lower panels: Ratio of the parallel (left) to perpendicular (right) mean square displacement between their three plasma modes. The subscripts i and j represent any two of Alfvén, fast, and slow modes. The simulations with the Larmor radius of Rg = 0.03Linj are based on the decomposed R1 data. The μ values for initially mirroring and scattering particles are set to μ0∈[0.1,0.2] and μ0∈[0.8,0.9], respectively.

In the text