The Diffusion and Scattering of Accelerating Particles in Compressible MHD Turbulence

We numerically study the diffusion and scattering of cosmic rays (CRs) together with their acceleration processes in the framework of the modern understanding of magnetohydrodynamic (MHD) turbulence. Based on the properties of compressible MHD turbulence obtained from observations and numerical experiments, we investigate the interaction of CRs with plasma modes. We find that (1) the gyroradius of particles exponentially increases with the acceleration timescale; (2) the momentum diffusion presents the power-law relationship with the gyroradius in the strong turbulence regime, and shows a plateau in the weak turbulence regime implying a stochastic acceleration process; (3) the spatial diffusion is dominated by the parallel diffusion in the sub-Alfvénic regime, while it is dominated by the perpendicular diffusion in the super-Alfvénic one; (4) as for the interaction of CRs with plasma modes, the particle acceleration is dominated by the fast mode in the high β case, while in the low β case, it is dominated by the fast and slow modes; and (5) in the presence of acceleration, magnetosonic modes still play a critical role in the diffusion and scattering processes of CRs, which is in good agreement with earlier theoretical predictions.

In general, the most classical acceleration mechanisms are considered as the second-and first-order Fermi processes (Fermi 1949;Bell 1978;Blandford & Ostriker 1978), and the (turbulent) magnetic reconnection (Sweet 1958;Parker 1957;Petschek 1964;LV99).Note that the second-order Fermi is also called stochastic acceleration process.This model, originally proposed by Fermi (1949), suggests that particles can statistically gain energy through collisions with interstellar clouds, which is similar to the reflection of particles due to magnetic mirror effects.Since MHD waves in the turbulence provide the motion of scattering centers, particles can be continuously scattered to advance (Melrose 1980).
As is well known, scattering is considered an essential process for CR acceleration.For instance, the scattering Gao & Zhang of CRs back into the shock is a necessary component of the first-order Fermi acceleration (see Longair 2011).At the same time, stochastic acceleration by turbulence is entirely based on the scattering process.Extensive numerical and analytical studies have been performed to understand the interactions of CRs with MHD turbulence, such as the scattering and diffusion processes (e.g., Yan & Lazarian 2002, 2004;Xu & Yan 2013;Lazarian & Xu 2021).When the pitch-angle1 approaches to 90 • , the scattering will vanish, and the mean free path becomes infinite (Fisk et al. 1974) since the particles are resonated by a very large wave-number k.This 90 • problem is one of the most famous concerns in quasi-linear theory's (QLT) predictions (Jokipii 1966).The root reason causing the 90 • problem is the assumption of unperturbed trajectories in QLT.In order to avoid this problem, Yan & Lazarian (2008) extends the QLT to nonlinear theory (termed NLT) by introducing finite resonance widths.
At present, there are a lot of simulation works focused on exploring the diffusion and scattering of particles in MHD turbulence.One part is associated with acceleration (e.g., Michalek et al. 1999;Beresnyak et al. 2011;Cohet & Marcowith 2016), while another with pure scattering and diffusion without acceleration (e.g., Xu & Yan 2013;Maiti et al. 2022;Hu et al. 2022).The works involving acceleration, are not comprehensive and have their limitations in application to a complex turbulent environment.Specifically, Michalek et al. (1999) only focused on a single turbulence regime, i.e., the cold plasma limit, where fast and slow (magnetosonic) waves degenerate to fast mode waves.Beresnyak et al. (2011) was based on the case of incompressible turbulence using MHD plus test particle simulations.In addition, Cohet & Marcowith (2016) focused on the sub-Alfvénic turbulence and explored the influence of driving ways on particle acceleration, where they found solenoidal forcing results do not match with QLT.Recently, Mertsch (2020) provided a more comprehensive overview of test particle simulations of CRs, which is based on QLT, its extensions, and the state-of-the-art in the test particle simulations mainly from the perspective of synthetic data simulation.
In the earlier studies, it was claimed that the Alfvén mode is inefficient for scattering and acceleration of CRs (Chandran 2000;Yan & Lazarian 2002) due to its anisotropic properties (Cho &Lazarian 2002, 2003, hereafter CL02 andCL03), while the fast mode is a major source of CRs scattering in the interstellar and intracluster media (Yan & Lazarian 2004;Brunetti & Lazarian 2007; see also Michalek et al. 1999 for the case of fast mode domination) due to its isotropic properties (CL02 and CL03).However, Zhang & Xiang (2021) demonstrated from the perspective of particle spectral distribution that the contribution of Alfvén mode to particle accelerations cannot be ignored, which even plays a dominant role in particle acceleration at the late stage of acceleration.In the case of incompressible turbulence, it has been claimed that the role of pseudo-Alfvén mode, i.e., slow mode, is crucial for the scattering of particles (Beresnyak et al. 2011;Xu & Lazarian 2020).
Furthermore, Demidem et al. (2020) claimed that the contributions of three modes to acceleration are comparable to within an order of magnitude in the relativistic MHD turbulence with Monte Carlo simulation of the test particle and synthetic data.It has been claimed that in the relativistic case, the properties of the three modes are similar to that of the non-relativistic one, except for the fast mode significantly coupling Alfvén mode (Takamoto & Lazarian 2016, 2017).Recently, Yuen et al. (2023) found that the pure Alfvén mode can be decomposed into an anomalous compressible component, which implies that the Alfvén mode may affect the particle transport.From the perspective of analytical research, in the framework of the modern understanding of MHD turbulence theory, Cho & Lazarian (2006) studied particle acceleration arising from fast and slow modes in both gaseous pressure-dominated (high β) and magnetic pressure-dominated (low β) cases.They predicted that the fast mode can accelerate particles more efficiently than the slow mode, and whether slow and fast modes dominate the acceleration of particles depends on the rate of spatial diffusion of particles.These interesting analytical predictions need to be tested and confirmed numerically.This motivates us to perform the current work by exploring various turbulence regimes.
One of the purposes of the current work is to study the diffusion and scattering behavior of particles being accelerated, and a second aim is to explore which plasma mode dominates the accelerated particle's transport in various turbulence regimes.We investigate how energetic particles are accelerated in MHD turbulence, and test which of the analytical predictions in Cho & Lazarian (2006) are consistent with our numerical results in a wide range of turbulence regimes that correspond to different astrophysical environments.Specifically, we want to know whether the fast mode dominates the acceleration, diffusion, and scattering process of particles exactly as the theoretical prediction, and whether the effect of the Alfvén mode is so weak as to be negligible.
This paper is organized as follows.In Section 2, we briefly introduce the theoretical description of MHD turbulence and diffusive properties of CRs.We perform test particle simulation and give our initial simulation set-up in Section 3. We present the numerical results, which contain the Mach numbers (Alfvénic Mach number M A , and sonic Mach number M s ) effect in four different turbulence regimes and the contributions of the three modes (wavelike isotropic fast mode, Alfvén, and slow Goldreich-Sridhar type modes) to particle acceleration and diffusion in Sections 4 and 5, respectively.Discussion and summary are provided in Sections 6 and 7, respectively.

MHD Turbulence Theory
The propagation of CRs is determined by their interactions with environmental turbulence.MHD turbulence theory has gone through a long period of development from pioneer works by Kolmogorov (1941;henceforth K41) and Iroshnikov & Kraichnan (1963;1965;henceforth IK65) to the modern MHD turbulence theory by Goldreich & Sridhar (1995;hereafter GS95) (see also Schekochihin 2020 for a recent review).By dimensional analysis, assuming that the energy is injected at a constant velocity v at large scales, K41 obtained the energy spectrum of E(k) ∝ k −5/3 for the pure fluid turbulence, called famous Kolmogorov spectrum.A few decades later, the IK65 spectrum, i.e., E(k) ∝ k −3/2 , was proposed for magnetized fluid turbulence.Its shortcoming is that it incompletely predicted an isotropic cascade of MHD turbulence.Later, GS95 suggested a scale-dependent anisotropy, i.e., k ∥ ≈ k 2/3 ⊥ L −1/3 for incompressible MHD turbulence, where L is the outer scale of the turbulence, and k ∥ and k ⊥ is the parallel and perpendicular components with regard to the local magnetic field of the wave vector k, respectively.
The modern theoretical understanding of MHD turbulence is based on the GS95 model.
In the framework of eddy motions, LV99 and Lazarian (2006) generalized incompressible MHD turbulence to compressible one in M A > 1 and M A < 1, respectively.In the case of M A < 1, with a sub-Alfvénic velocity driving turbulence at injection scale L inj , the weak turbulence cascade spans the range from L inj to l tr , where l tr = L inj M 2 A is the transition scale at which we have turbulence velocity v l equal to Alfvén speed V A , while the strong turbulence cascade is in the range of [l dis , l tr ], where l dis is the dissipation scale.And in this range, we have the relationship of associated with the parallel scale l ∥ and transversal one l ⊥ , which suggests that eddies are stretched along the local magnetic field.Note that the original GS95 relation of l ∥ ∝ l 2/3 ⊥ can be recovered by setting M A = 1 in Equation (1).
As for M A > 1, i.e., the super-Alfvénic turbulence, the cascade starting from the injection scale L inj has a hydrodynamic Kolmogorov property, because of the marginal influence of a weak magnetic field on the cascade process.When the cascade is decreased to l tr = L inj M −3 A (Lazarian 2006), it enters a regime of strong turbulence, from which to the dissipation scale it again exhibits the GS95 anisotropy characteristics, with the scale relation of l ∥ ≈ L 1/3 inj l 2/3 ⊥ M −1 A and the velocity-scale relation of v l = V L (l ⊥ /L inj ) 1/3 , where V L is the turbulence injection velocity at the scale L inj .
Furthermore, the analytical solution of the MHD dispersion equations can give three solutions (e.g., Beresnyak 2019) that correspond to Alfvén, fast, and slow modes as numerically confirmed by CL02.Among the three modes, Alfvén mode is incompressible, while slow and fast modes are compressible and are called magnetoacoustic waves/modes.
The latter two have scale-dependent anisotropic properties of l ∥ ∝ l 2/3 ⊥ in the local frame of the magnetic field (LV99; Cho & Vishniac 2000;Maron & Goldreich 2001).Slow mode, passively mixed by Alfvén mode (Lithwick & Goldreich 2001), has the same anisotropic scaling as Alfvén mode.Therefore, Alfvén and slow modes present the K41 spectrum of

Physical Quantities That Characterize CR Propagation
In general, the propagation and acceleration of CRs can be described by the Fokker-Planck equation coefficients (see the book from Schlickeiser 2002).Here, we mainly focus on the coefficients relevant to our following simulations, such as the pitch-angle, spatial, and momentum diffusion coefficients.From a theoretical point of view, the diffusion of CRs results from the resonant (gyroresonance) and nonresonant (transit time damping, i.e., TTD) interaction of CRs with MHD turbulence.However, from the perspective of numerical simulation, it is currently difficult to distinguish the contribution of individual components.
As for the pitch-angle scattering, the pitch-angle diffusion coefficient is defined by where the square brackets indicate an average over the ensemble of particles at the integration time t.The pitch angle θ is the angle between the particle velocity vector and the local magnetic field direction, whose cosine value is µ = cosθ at the time t, while µ 0 is at the initial moment t 0 .
Here, pitch angle θ range is from 0 to 90 • , which corresponds to both µ and µ 0 from 1 to 0. In the QLT, the mirror resonance has a sharp peak at 90 • due to discrete Landau resonant condition, resulting in the disappearance of the diffusion coefficient (Goldstein 1976;Felice & Kulsrud 2001).This is the infamous 90-degree problem.To avoid the 90 • problem, QLT was extended to NLT by taking the magnetic mirroring effect into account on large scales (Yan & Lazarian 2008).
In addition, when introducing the bouncing of CRs, this problem of quasilinear gyroresonant can be alleviated (see Lazarian & Xu 2021 for the details).
Using the pitch-angle diffusion coefficient, we can determine the parallel mean free path of the particles by substituting D µµ into (Earl 1974) where u is the velocity of particles.Alternatively, the mean free path λ ∥ of particles can be calculated by Here, D ∥ is a parallel diffusion coefficient (Giacalone & Jokipii 1999) where ( x − x0 ) is the spatial separation in the local magnetic field directions.Numerically, it is not difficult to test the fact that Equations ( 4) and ( 3) give similar results (see also Maiti et al. 2022 for their testing).Similarly, we can also obtain the perpendicular diffusion coefficient where (ỹ − ỹ0 ) represents the spatial separation perpendicular to the local magnetic fields.
Assuming that the particles move in terms of a random walk in momentum space, we can determine the momentum diffusion coefficient averaged over the ensemble of particles, where ∆p is the amount of change in particle momentum in the time interval ∆t.In this work, we use Equation ( 7) to characterize the particle diffusion behavior in the acceleration processes.

Simulation of MHD Turbulence
The third-order-accurate hybrid, essentially non-oscillatory code is adopted to solve the ideal MHD equations describing MHD turbulence as follows: where p g = c 2 s ρ is the gas pressure, c s and ρ represent the sonic speed and density, respectively, t is the time of fluid evolution, J = ∇ × B is the current density, and f is a random solenoidal driving force on large scale (small wavenumber k ≃ 2.5 in wavenumber space) in our simulation.The other parameters have their usual meaning.Our simulation is performed in a periodic box at the length of L = 2π, setting a non-zero mean magnetic field strength along the x-axis direction.When statistical steady-state is reached, we output primitive physical quantities such as the magnetic field, velocity, and density.To characterize each simulation, we calculate Alfvénic number by √ ρ is the Alfvén speed.The resulting values are listed in Table 1, where each model corresponds to different turbulence regimes.
Numerically, the compressible MHD turbulence was first decomposed into Alfvén, slow and fast modes by a Fourier transformation (CL02; CL03).The limitation of this method is that it only applies to the global reference frame and can only deal with the problems of M A < 1.Later, Kowal & Lazarian (2010) improved this separation technique by introducing a discrete wavelet transformation before the Fourier separation, extending it to the M A > 1 case.The displacement vectors of the slow, fast, and Alfvén modes are defined by their unit vectors where ϕ is an angle between k and B 0 .After numerically determining the ζs , ζf and ζA , we can project the total magnetic field and velocity into these three unit vector directions and obtain information about the magnetic and velocity fields of each mode.To maintain consistency in the operation of all data, we use the wavelet decomposition in this work.

Method of Test Particle
In the electromagnetic field, the kinetic equation that the motion of a charged particle satisfies is where γ ≡ 1/ 1 − u2 /c 2 is the Lorentz factor of the particle, m and q are the mass and charge of the particle.In this work, we take into consideration the acceleration process resulting from magnetic field B and its inductive field E = −v×B, without considering the resistivity of the electric field.Therefore, Equation ( 15) can be rewritten as By integrating this equation, we can trace particle's trajectories using the classic 4th order explicit Runge-Kutta (RK4) method with adaptive time step. 2 We can recover the local values of the plasma velocity v and magnetic field B at each step of the integration by using linear interpolation with the discontinuity detector.
When analyzing the numerical output, we calculate the direction of the local magnetic fields by (Cho & Vishniac 2000) And then we can define D ∥ as the parallel direction of x = B l /|B l | and D ⊥ as its perpendicular direction.

NUMERICAL RESULTS: THE INFLUENCE OF TURBULENCE PROPERTIES ON PARTICLE ENERGIZATION
In this section, we provide numerical results of particles' energization processes related to different turbulence regimes, based on data cubes from simulations of MHD turbulence listed in Table 1.Specifically, our results include the trajectory and gyroradius of particles, their momentum, spatial and pitch-angle diffusion in four turbulence regimes such as sub-Alfvénic & subsonic, sub-Alfvénic & supersonic, super-Alfvénic & subsonic, and super-Alfvénic & supersonic.

Trajectory of Particles
At the beginning of the simulation, we instantly inject 10 4 test particles with a Maxwell-type spectrum randomly distributed into the whole box space of the simulation.Throughout this paper, we set the temperature of test particles to be 10 6 K, the plasma velocity (i.e., the fluid velocity) to be 5 × 10 7 m/s = 0.1667 c, and magnetic field strength to be 5 µGs.
Figure 1 plots the trajectories of the 1000th, 5000th, and 7000th test particles.As shown, the particles experience diffusive motion in the simulation space.We see that the particles cross a large spatial scale in the x-axis direction, which is because, in our simulation of MHD turbulence, the mean magnetic field is set along this direction.Compared with panels (a) and (b), we see also that the particles in panels (c) and (d) span a larger space in three coordinate axis directions.Given that the setting of the magnetic field strength of 5.0 µGs, we have the dimensionalized mean magnetic field strength of B 0 = 5.0 µGs for sub-Alfvénic turbulence (panels a and b) and of B 0 = 0.5 µGs for super-Alfvénic turbulence (panels c and d).According to the formula of gyro motion of a charged particle, where p ⊥ = γmv ⊥ is the perpendicular momentum of a particle at the perpendicular velocity v ⊥ with respect to the local magnetic field, B, it is not difficult to understand particles' motion behavior.For a fixed p ⊥ , the larger the B value is, the smaller the R g .As a result, particles have a smaller gyroradius, resulting in a smaller spatial extension (upper panels).Another interesting point is that in the case of super-Alfvénic & subsonic turbulence (lower-left panel), particles have the largest gyroradii at the final acceleration time (see also the lower-left panel of Figure 2), gaining the maximum kinetic energy of E max = 715.73m p c 2 .This implies that a strong mean magnetic field does not necessarily result in the maximum possible energy.Except for tracking the trajectory of a single particle, we also check particle's drift motion by observing the evolution of particle velocities and displacement over time (not shown here).We find that the perpendicular velocity with respect to the local mean magnetic field is always greater than the parallel one.As for the time evolution of the displacement, we find that the perpendicular displacement is much larger than the parallel one in the early stage of the evolution, which indicates that magnetic field gradient drift may dominate the acceleration.However, in the late stage, the larger parallel displacement indicates the domination of turbulence cascade interaction, which is consistent with the results from diffusion analysis (see Section 4.3).

Momentum Diffusion
Figure 2 plots the evolution of the gyroradius over time for four turbulence regimes, where the red solid lines represent the averaged gyroradius of all particles and the cyan dashed lines represent the gyroradius distribution of a single particle.
As seen, the evolution of the gyroradius over time presents an overall consistency in the different turbulence regimes.Based on the length scales, we divide each evolution of R g into three main stages.In the first stage (R g < L/512), the gyroradius R g of particles slowly increases with the time after a plateau.This indicates that at the start of the simulation, the acceleration of the particles is inefficient, which implies that the particles cannot interact with turbulence waves (scattering centers) efficiently.The acceleration phenomenon of R g less than one grid size (L/512) may arise from the non-resonant interaction between the particles and fluid, the gradient and curvature drifts of magnetic fields.
In the second stage (L/512 < R g < L tr ), we see that R g shows a power-law relationship of t 4/3 for the subsonic turbulence, while for the supersonic turbulence, R g shows the power law of t 5/3 and t 1.0 in the sub-Alfvénic and super-Alfvénic regimes, respectively.In this strong turbulence regime we are mostly interested in, the acceleration efficiency of particle acceleration has been significantly improved, with a distribution ranging from R g ∝ t 1.0 to R g ∝ t5/3 .It arises from quasi-resonant interactions of particles with different scale eddies, that is, the interaction between fluid and magnetic fields causes the induced electric field to accelerate particles.In the third stage (R g > L tr ), the time evolution of R g for all four regimes shows a power law of R g ∝ t 2/3 in this weak turbulence regime.This demonstrates that the efficiency of particle interaction with turbulence is decreased.
Figure 3 shows the evolution of the momentum diffusion coefficient over time (left column) and gyroradius (right column), including the total momentum diffusion D pp (upper row), its perpendicular component D pp,⊥ with respect to the local magnetic field (middle row), and the ratio of parallel to perpendicular components D pp,∥ /D pp,⊥ (lower row). 3As shown in panel (a), before 1/Ω 0 , corresponding to the first stage of R g < L/512 mentioned above, the total D pp over time presents a power-law relationship of D pp ∝ t4/5 .At t ≃ 1/Ω 0 , the evolution of D pp over time presents a trough (seen also in Demidem et al. (2020)) except for the super-Alfvénic and supersonic one (orange-dotted line), which indicates that the turbulence interaction tends to suppress the particle diffusion in the dissipation region.When entering the strong turbulence regime, the increasing rate of D pp for sub-Alfvénic turbulence gets higher with a steeper distribution.In these two stages, the particle undergoes the superdiffusion in the momentum space 4 .When the gyroradius reaches the transition scale, i.e., enters the weak turbulence regime, the evolution of D pp follows a plateau phase, which is a significant feature of the second-order Fermi process (Pezzi et al. 2022;Liang et al. 2023) and means that the particle experiences the normal diffusion.
To explore the relationship between D pp and the gyroradius R g , we divide R g from 10 −3 /L to 10 2 /L into 50 bins in the logarithmic bin, then pick out the particles in each bin to calculate their diffusion coefficients. 5In Gao & Zhang the box length range from R g = L/512 to L that we are interested in, we plot in panel (d) the momentum diffusion coefficient as a function of R g .As shown, it approximates to D pp ∝ R g 3/4 in the strong turbulence regime for sub-Alfvénic turbulence, while D pp ∝ R g 2/5 for super-Alfvénic turbulence.This demonstrates that momentum diffusion in sub-Alfvénic turbulence is faster than that in super-Alfvénic turbulence.This should be the fact that in the momentum space, the strong magnetic fields improve particle diffusion.For both cases of sub-and super-Alfvénic turbulence, D pp shows a plateau after R g > L tr , i.e., after 10/Ω 0 ; see also panel (a).Furthermore, we also explore the evolution of D pp with momentum p, which shows a similar power-law relationship with that of gyroradius, i.e.D pp ∝ p 3/4 for sub-Alfvénic turbulence, and D pp ∝ p 2/5 for super-Alfvénic one.Note that Cho & Lazarian (2006) predicted D pp ∼ (∆p) 2 /∆t ∼ p 2 (∇•v l ) 2 ∆t, from which we can reasonably speculate ∇•v l ∝ p 5/8 /∆t and ∝ p 4/5 /∆t.
In the middle panels of Figure 3, we plot the perpendicular component D pp,⊥ over the integrated time and gyroradius (by binning the particles into different values of gyroradius).Its overall behavior is similar to D pp , except for some small differences.In particular, we see a good power-law relation of D pp,⊥ ∝ t 4/5 before 1/Ω 0 and D pp,⊥ ∝ R g 6/5 in the strong turbulence regimes for four turbulence scenarios.Moreover, the ratio of D pp,∥ and D pp,⊥ is shown in the lower panels, from which we can see that at the early stage of evolution D pp,⊥ is even one order of magnitude larger than D pp,∥ for the sub-Alfvénic & super-Alfvénic turbulence regime.This implies that at small R g the perpendicular gradient drift of the magnetic field dominates particle diffusion processes by v grad = γmv 2 ⊥ (B × ∇B)/2q 2 B 2 .For the time evolution of D pp (panel c), after t ∼ 1/Ω 0 , the diffusion is dominated by the parallel momentum.For its gyroradius evolution (panel f), in the case of super-Alfvénic turbulence, when R g > L/512, D pp,∥ /D pp,⊥ grows slowly and enters the stage dominated by parallel momentum, and in the sub-Alfvénic case, when R g > L tr the parallel momentum dominates the momentum diffusion.This means that the diffusion is dominated by the parallel momentum at large R g .

Spatial Diffusion and Pitch-Angle Scattering
We bin particles into different values of gyroradius, then present the spatial diffusion coefficient as a function of the gyroradius in Figure 4, including the parallel diffusion D ∥ (panel a) and the ratio of parallel and perpendicular components (panel b) for four turbulence regimes.In panel (a), D ∥ for sub-Alfvénic turbulence shows a power-law relationship of ∝ R 1/4 g and a plateau before and after R g ∼ 1/L, respectively.This indicates that the spatial diffusion of particles first increases in the strong turbulence regime, i.e., scales smaller than L tr , and then enters a plateau stage in the weak turbulence regime, i.e., scales greater than L tr .In the case of super-Alfvénic turbulence, D ∥ shows a slower slope at scales smaller than L tr , which implies that strong magnetic fields enhance the parallel spatial diffusion.
To gain insight into the detailed information of the spatial diffusion of the accelerated particles, we further show the ratio of parallel and perpendicular diffusion coefficients D ∥ /D ⊥ in panel (b).In the case of M A < 1, the ratio is larger than 1 all the time, which means that D ∥ dominates the spatial diffusion in this case.As for the case of M A > 1, the distribution of the ratio can be divided into two stages: first, it is larger than 1 and slowly decreases, which means that the spatial diffusion of CRs is dominated by D ∥ ; second, after R g ∼ 0.3/L, the ratio is less than 1 and enters into the stage dominated by D ⊥ .This indicates that in the strongly turbulent cascade process, super-Alfvénic turbulence tends to inhibit the parallel diffusion of the accelerated particles.Similarly, by analyzing the mean free path of the particles (see Equation 4), we also come to similar results.
Next, with the aim of studying the pitch-angle scattering behavior of the accelerated particles, we bin particles into different values of R g and show the evolution of D µµ over R g in Figure 5.With increasing R g , D µµ can be roughly divided into two stages.Firstly, differences between these four turbulence regimes begin to emerge during the plateau stage of R g vs. t (see Figure 2).D µµ of these four turbulence regimes presents a similar evolution with R g , which decays with a power-law relationship of D µµ ∝ R g −1/5 , except for a slightly slower decay of the sub-Alfvénic and supersonic case (green triangle-down line).Secondly, when R g of particles reaches the transition scale L tr , where the momentum and spatial diffusion reach a plateau, the decay rate of D µµ slows down and gradually enters a plateau for these four regimes.Note that there are opposite evolutionary trends between D µµ and D pp over R g .This suggests that when a particle is being accelerated, its pitch-angle scattering is suppressed.

Trajectory of Particles
Based on data cubes A1 and A5 listed in Table 1 from simulations of MHD turbulence, we provide numerical results of CRs interacting with decomposed plasma modes.Figure 6 plots the trajectories of selected three particles from the pre-decomposed MHD turbulence (panel a) and its post-decomposed three modes: the Alfvén (panel b), fast (panel c), and slow (panel d) modes.The pre-decomposed MHD turbulence shows a similar diffusive motion as Figure 1.As for the post-decomposed plasma modes, the trajectories for the fast mode are similar to that of pre-decomposed turbulence, and trajectories for the Alfvén mode are extended along the x-axis direction, which is consistent with our expectation due to its intrinsically polarized feature.Although the slow mode has similar trajectories with the Alfvén one, which should be from their same scale-dependent anisotropies (Cho & Lazarian 2002), it experiences more extended spatial lengths.This phenomenon implies that the particle interaction with turbulence for fast mode is more effective, and the particle acceleration rate should be higher in the supersonic regime.In addition, we see that the spatial scales in the x-axis direction for slow and fast modes are slightly larger than that for Alfvén mode, which may mean that magnetosonic modes have a slightly higher diffusion efficiency than Alfvén mode.These inferences will be further confirmed below.In Figure 7, we show the gyroradius of particles from the three modes and their ratios in the high and low β regimes. 6s shown in the left column for the high β regime, we can see that the acceleration processes can be divided into three different stages by the length scales.In the first stage (R g < L/512), the gyroradii of these three modes present a plateau, and then slowly increase over time to exceed L/512.In the second stage (L/512 < R g < L tr ), i.e., the strong turbulence regime, R g of fast and Alfvén modes is very similar to each other presenting a power-law relationship of t 4/5 , and larger than R g,S at the same simulation time.This may be related to the fraction of magnetic energies listed in Table 1, that is, the larger the B value, the R g is smaller (see Equation ( 18)).In the third stage (R g > L tr ), i.e., the weak turbulence regime, the fast mode remains the same power law as that in the second stage, while Alfvén and slow modes show a shallower power law of R g ∝ t 1/2 .Interestingly, the acceleration behavior of particles caused by three plasma modes is different from that caused by overall MHD turbulence before being decomposed (see Figure 2).

Momentum Diffusion
To explore their ratios, we highlight the timescale zone in panel (b) that the gyroradii of three modes reach the L/512 and L tr scales, by filling in the area between two adjacent vertical lines at each scale for convenience.We can more clearly see their differences: during the plateau stage (i.e., the left region of the pink vertical bandwidth), R g,A ≃ R g,F ≃ R g,S ; after that (i.e., the region between the pink and cornflower-blue vertical bandwidth), R g,S is almost always smaller than R g,A and R g,F ; as for the weak turbulence regime (i.e., the right region of the violet vertical bandwidth), it is opposite that R g,F is the largest one then followed by R g,A .
In the case of the low β regime (right column), there is a distinctly different phenomenon for the slow mode before about 10 0.5 /Ω 0 (during the strong turbulence regime).Slow wave starts to accelerate before the plateau stage of fast and Alfvén modes.After 10 0.5 /Ω 0 , R g of the slow mode grows slowly and presents a power-law relationship of R g ∝ t 2/5 , compared with the high β scenario.As seen, in the first two stages (the left region of the corn-flower-blue vertical bandwidth), the interaction caused by the slow mode dominates the acceleration process, while in the third stage (the right region of the corn-flower-blue vertical bandwidth), the acceleration efficiency of the fast mode is higher than that of the other two.In this comparison study, since Models A1 and A5 listed in Table 1 have very close values of M A = 0.65 and 0.5, respectively, the difference in β is mainly caused by the difference in M s (0.48 for A1 and 9.92 for A5) by the Figure 6.The 3D trajectories of three test particles selected, arising from data A5 (panel a) listed in Table 1 and its decomposed modes including the Alfvén (panel b), fast (panel c), and slow (panel d).The star and circle markers refer to their initial and final locations in the 3D space, respectively.L box = 512 is the box scale.
relationship of β = 2M 2 A /M 2 s .Therefore, we find that the shock wave occurring in the supersonic case influences the particle's acceleration.
To explore the effect of the plasma modes on the diffusion in momentum space, we show in Figure 8 the ratio of momentum diffusion coefficients between two modes as a function of the time, for the total momentum (top row), its parallel D pp,∥ (middle row), and perpendicular D pp,⊥ (bottom row) components.In the high β regime (left column), for the total momentum diffusion coefficient (panel a), the contribution of magnetosonic modes (D pp,S and D pp,F ) is greater than Alfvén mode (D pp,A ) during the whole evolution.As for D pp,∥ (panel b) and D pp,⊥ (panel c), it is basically the same as the total momentum diffusion coefficient, though D pp,F /D pp,S is close to one after about 1/Ω 0 , corresponding to the last two stages of the time evolution of gyroradius.This means that in the high β case, the fast mode is as crucial as the slow mode for the parallel and perpendicular components of momentum diffusion in the later stages of the particle acceleration.
In the low β regime (right column), the contribution of the fast mode is almost always the largest one, followed by Alfvén mode, and slow mode is the smallest, including the total momentum diffusion coefficient D pp , as well as its parallel and perpendicular components (D pp,∥ and D pp,⊥ ).In detail, the order of magnitude of the ratio of any two modes in low β is slightly higher than that in high β.Especially for D pp,∥ , the fast mode is almost several orders of magnitude higher than the slow mode in the first stage.
In short, our one finding is that the fast mode in the high β case dominates the particle acceleration, while in the low β case, the fast and slow modes dominate the acceleration.Here, the dominance of the fast mode is in agreement with the earlier theoretical predictions (Yan Another finding is that our numerical results confirm the dominance of magnetosonic modes in the acceleration and momentum diffusion of CRs, which is consistent with the previous studies (Yan & Lazarian 2004;Lynn et al. 2014;Zhang et al. 2020).

Spatial Diffusion and Pitch-Angle Scattering
To explore the effect of the plasma modes on particle spatial diffusion, we present the ratio of parallel D ∥ (upper row) and perpendicular diffusion coefficient D ⊥ (lower row) between two modes in Figure 9.As is shown in the left column (high β), the slow mode dominates the parallel diffusion D ∥ (panel a) during the whole simulation time, except for the time range after ∼ 10/Ω 0 (D ∥,A ≈ D ∥,F ≈ D ∥,S ).From panel (b), we see that the dominance of D ⊥ among three modes can be roughly divided into two stages: (1) slow mode dominates perpendicular diffusion followed by D ⊥,F > D ⊥,A before ∼ 10/Ω 0 ; (2) after 10/Ω 0 , magnetosonic modes dominate perpendicular diffusion with the relation of D ⊥,F ≈ D ⊥,S ≈ D ⊥,A .It should be noted that the Alfvén mode plays a sub-dominant role almost in the whole evolution.
In the case of the low β (right column), the ratio of parallel and perpendicular diffusion coefficients is almost similar.Before ∼ 0.1/Ω 0 , both parallel D ∥ and perpendicular D ⊥ keep a dominant order, that is, fast mode is the largest followed by Alfvén mode, and slow mode is the smallest.From ∼ 0.1/Ω 0 to ∼ 10 2 /Ω 0 , the evolution of both D ∥ and D ⊥ in this situation is similar to the case of D ∥ in high β, where magnetosonic modes dominate the diffusion.After ∼ 10 2 /Ω 0 , their ratios are all close to 1, which means that the three modes play a comparable role.As a result, magnetosonic modes play a dominant role in the spatial diffusion of particles in the strong turbulence regime.
D µµ for three modes in Figure 10.As is shown in panel (a), i.e., the case of high β, the dominant relation is D µµ,S > D µµ,F ∼ D µµ,A before 1/Ω 0 and then it turns to the stage that three modes approximately keep the similar scattering level.Differently, for the low β regime (panel b), the fast mode is the largest, and Alfvén mode is secondary before ∼ 0.25/Ω 0 , corresponding to R g reaching the scale of L/512 Figure 8.The ratio of momentum diffusion coefficients between two modes as a function of the evolution time arising from total momentum (top row), and its parallel (middle row) and perpendicular (bottom row) components, based on A1 (left column) and A5 (right column) as listed in Table 1.The light curves with many fluctuations are original distributions, while the dark curves are a robust local weighted regression fitting.The other descriptions are the same as those of Figure 7.
for fast and Alfvén modes (see also Figure 7).After that, three ratios are approximately equal to 1 (after presenting a trough) consistent with the high β scenario (see panel b).We would like to stress that when plotting, we constrain a narrow range of vertical coordinates to observe their differences.Therefore, using pitch-angle diffusion versus the evolution time, it is difficult to distinguish their differences from the scattering processes of the accelerated particles in the range of the box size.

DISCUSSION
In this work, we numerically explore the energization processes of CRs in compressible MHD turbulence regimes.We first focus on the influence of different turbulence properties on the CR energization and then on the interaction of CRs with Alfvén, slow, and fast modes.Specifically, we elaborate on the acceleration, diffusion, and scattering behavior of 10 4 particles injected instantly by considering their evolution with the time and gyroradius.
Our current studies are built on a modern understanding of MHD turbulence theory, focusing on the properties of turbulence and the influence of plasma modes on particle acceleration.The results demonstrate that the diffusion coefficients of accelerated particles experience a universally exponential increase with the simulation time and gyroradius, similar to the phase of exponential growth seen in Demidem et al. (2020) for the relativistic MHD turbulence.Differently, Demidem et al. (2020) used the Monte Carlo simulation, one of the merits of which is that it is not affected by any choice of the distribution function.However, the limitation of this method is that its phenomenological description of acceleration processes When exploring the interplay of CRs and individual plasma modes, we find that the contribution to acceleration arising from the fast mode is the highest one in the high β case, while the fast and slow modes dominate the acceleration in the low β one, which is in good agreement with Cho & Lazarian (2006) for the particle acceleration under strong and weak MHD turbulence.They theoretically derived the momentum diffusion coefficients of fast and slow modes and established that the acceleration efficiency of fast mode is more efficient than that of slow mode similar to the results of Yan & Lazarian (2002) and Chandran (2003).The earlier studies have claimed that the scattering of charged particles by Alfvénic turbulence is negligible (Yan & Lazarian 2002, 2004;Lynn et al. 2014), and magnetosonic modes are dominant processes for the transport and acceleration of CRs.However, we found that the effect of Alfvén mode cannot be ignored, playing a sub-dominant role.We conjecture that the non-negligible impact of the Alfvén mode on particle acceleration and diffusion may be from the compressible component of the Alfvén waves.In this regard, by projecting the local Alfven component of the turbulent variables into linear combinations of Alfvén and non-Alfvén components, Yuen et al. (2023) proposed an Alfvén leakage effect in the frame of the local magnetic field.
It is generally accepted that the mechanism by which turbulence causes particle acceleration is a stochastic acceleration process.Our current work does not describe what the acceleration mechanism is.As mentioned above, this paper focused on the influence of both turbulence properties and plasma modes on the transport of CRs.As well-known, the interaction of MHD turbulence would result in the gyroresonance and TTD of particles.For the former, its criterion is based on a comparison between wave frequency and particle Larmor frequency.At small scales, a particle can preserve its adiabatic invariant mv 2 ⊥ 2B 0 = const, due to the Larmor radius being smaller than the variation scale of the magnetic field.However, particle motions would violate the adiabatic invariant conservation at large scales (Chandran 2000;Yan & Lazarian 2003;Lazarian & Xu 2021).For the latter, it is essentially a Cerenkov-type interaction allowing particles to interact with large-scale turbulence.In the framework of the modern understanding of MHD turbulence, separating out the contributions of TTD and gyroresonance will be performed in future work.
When visualizing the time evolution of the gyroradius of particles (and also including the interaction of CRs with decomposed plasma modes), we use the weak-strong transition scale from theoretical predictions in LV99 and Lazarian (2006) to explain our numerical findings.The theoretical prediction of weak-strong transition scales has been confirmed by simulations (Verdini & Grappin 2012;Meyrand et al. 2016;Makwana & Yan 2020) and observations (Sioulas et al. 2023;Zhao et al. 2023).As an example, in the perspective of simulation, Meyrand et al. (2016) presented direct evidence of such a weak-strong transition, using a high-resolution three-dimensional direct numerical simulation of incompressible MHD turbulence.From the perspective of observation, Zhao et al. (2023) reported the first observational evidence for the Alfvénic weak-strong transition in MHD turbulence in Earth's magnetosheath using data from the Cluster spacecraft.
In order to verify the reliability of our results, we also perform in this paper a comparison study using the numerical resolutions of 256 3 and 792 3 .Our studies show that the difference in the numerical resolution does not affect the simulation results provided in the current work.Therefore, the resolution of 512 adopted in this paper is sufficient for our current goal.This work is only a first step towards understanding a more complex interplay between particle acceleration and plasma modes.

SUMMARY
This paper studied the interaction of CRs with compressible MHD turbulence together with their acceleration processes based on the modern understanding of MHD turbulence theory.With different MHD turbulence regimes that may happen in a realistic astrophysical environment, we focused on particle acceleration, diffusion, and scattering processes using the test particle simulation.
1. We find that the gyroradius of particles exponentially increases with the simulation time by R g ∝ t φ .In the strong turbulence regime, the index φ ∈ [4/3, 5/3] for sub-Alfvénic turbulence is steeper than φ ∈ [1, 4/3] for super-Alfvénic one, while in the weak turbulence regime, φ is approximately equal to 2/3 for four turbulence cases explored.
2. In the strong turbulence range, the particle undergoes the superdiffusion in the momentum space, with the relationships of D pp ∝ R 3/4 g for sub-Alfvénic turbulence and D pp ∝ R 2/5 g for super-Alfvénic one.The momentum in the direction parallel to the local magnetic field dominates the diffusion process at large R g .In the weak turbulence regime, the momentum diffusion shows a plateau implying a stochastic acceleration process and meaning that the particle experiences the normal diffusion in the momentum space.
3. With R g in the range of box size, the parallel diffusion dominates the spatial diffusion of particles in the case of the sub-Alfvénic turbulence regime, while in the weak turbulence regime the perpendicular diffusion is slightly faster than the parallel one in the case of super-Alfvénic one.The pitch-angle diffusion decays along with the increasing gyroradius, with the diffusion rate in the weak turbulence slower than that in the strong turbulence.
4. As for the interaction of CRs with individual plasma modes, the property of particle acceleration, diffusion, and scattering is distinct from that of pre-decomposed MHD turbulence.

The particle acceleration is
Gao & Zhang dominated by the fast mode in the high β case, while in the low β case, it is dominated by the fast and slow modes.Moreover, magnetosonic (fast and slow) modes are also the main contributor to the momentum diffusion of CRs.
5. The spatial diffusion of particles is dominated by the slow mode for both the cases of high and low β in the strong turbulence regime, while in the weak turbulence regime, three plasma modes play a comparable role.In particular, the spatial diffusion from the Alfvén mode cannot be ignored.
We would like to thank the anonymous referee for constructive comments that have significantly improved our manuscript.We thank Alex Lazarian for reading throughout the manuscript, and Siyao Xu for helpful discussions on the particle acceleration and scattering in MHD turbulence.
The authors thank the support from the National Natural Science Foundation of China (grant No. 11973035).J.F.Z. also thank the Hunan Natural Science Foundation for Distinguished Young Scholars (No. 2023JJ10039), and the China Scholarship Council for overseas research fund.

Figure 1 .
Figure 1.The 3D trajectories of three test particles selected, arising from sub-Alfvénic & subsonic (panel a), sub-Alfvénic & supersonic (panel b), super-Alfvénic & subsonic (panel c), and super-Alfvénic & supersonic (panel d) turbulence regimes.The star and circle markers refer to their initial and final locations in the 3D space, respectively.L box = 512 is the box scale.

Figure 2 .
Figure 2. The time evolution of gyroradius of a single particle is compared with the behavior of an ensemble of 10 4 test particles interacting with four different turbulence scenarios.The horizontal dashed, dotted and dash-dotted lines show the grid, transition and box scales, respectively.L and Ω 0 are the box length and initial gyrofrequency, respectively.

Figure 3 .
Figure 3.The momentum diffusion coefficient as a function of the evolution time (left column) and the gyroradius of particles accelerated (right column).The total momentum diffusion coefficient, its perpendicular component, and the ratio of parallel to perpendicular component are plotted in the upper, middle, and lower rows, respectively.The vertical dashed and dash-dotted lines plotted in the right column present the grid and box sizes, respectively, while the horizontal dash-dotted line shows D pp,∥ /D pp,⊥ = 1.L, Ω 0 , and P 0 are the box length, initial gyrofrequency, and momentum, respectively.

Figure 4 .
Figure 4.The evolution of parallel diffusion coefficient (panel a), D ∥ , and the ratio of D ∥ to perpendicular diffusion coefficient D ⊥ (panel b) with the gyroradius.The other descriptions are the same as those of Figure 3.

Figure 5 .
Figure 5.The pitch-angle diffusion coefficient, D µµ , as a function of gyroradius at four turbulence regimes explored.The other descriptions are the same as those of Figure 3.

Figure 7 .
Figure 7.The gyroradius of the accelerated particles as a function of the evolution time for decomposed three plasma modes (upper row) and the ratio of three modes (lower row).The left and right columns are based on A1 and A5 listed in Table 1, respectively.The subscript i represents Alfvén and fast modes, while the subscript j slow and fast modes.The vertical dashed and dotted lines represent the time that the gyroradius reaches the grid and transition scales, respectively, where the colors yellow-green, orange, and tomato present the Alfvén, fast, and slow modes, respectively.L and Ω 0 are the box length and initial gyrofrequency, respectively.& Lazarian 2002; Chandran 2003; Cho & Lazarian 2006).Another finding is that our numerical results confirm the dominance of magnetosonic modes in the acceleration and momentum diffusion of CRs, which is consistent with the previous studies(Yan & Lazarian 2004;Lynn et al. 2014;Zhang et al. 2020).

Figure 9 .
Figure 9.The ratio of parallel (upper row) and perpendicular (lower) diffusion coefficients for three modes.The other descriptions are the same as those of Figure 8. cannot trace the microphysical details related to plasma physics.When exploring the interplay of CRs and individual plasma modes, we find that the contribution to acceleration arising from the fast mode is the highest one in the high β case, while the fast and slow modes dominate the acceleration in the low β one, which is in good agreement withCho & Lazarian (2006)  for the particle acceleration under strong and weak MHD turbulence.They theoretically derived the momentum diffusion coefficients of fast and slow modes and established that the acceleration efficiency of fast mode is more efficient than that of slow mode similar to the results ofYan & Lazarian (2002) andChandran (2003).The earlier studies have claimed that the scattering of charged particles by Alfvénic turbulence is negligible(Yan & Lazarian 2002, 2004;Lynn et al. 2014), and magnetosonic modes are dominant processes for the transport and acceleration of CRs.However, we found that the effect of Alfvén mode cannot be ignored, playing a sub-dominant role.We conjecture that the non-negligible impact of the Alfvén mode on particle acceleration and diffusion may be from the compressible component of the Alfvén waves.In this regard, by projecting the local

Figure 10 .
Figure 10.The ratio of pitch-angle diffusion coefficients for three modes.The other descriptions are the same as those of Figure 8.

Table 1 .
Data cubes with numerical resolution of 512 3 for different turbulence regimes.M A and M s are the Alfvénic and sonic Mach numbers, respectively, and β = p gas /p mag = 2M 2