Energization of Charged Test Particles in Magnetohydrodynamic Fields: Waves versus Turbulence Picture

Direct numerical simulations of three-dimensional compressible magnetohydrodynamic (MHD) turbulence have been performed in order to study the relation between wave modes and coherent structures and the consequent energization of test particles. Moreover, the question of which is the main mechanism of this particle energization is rigorously discussed. In particular, using the same initial conditions, we analyzed the nonlinear and linear evolution of a turbulent state along with the case of randomized phases. Then, the behaviors of the linear and nonlinear simulations were compared through the study of the time evolution of particle kinetic energy and preferential concentration. Also, spatiotemporal spectra were used to identify the presence of wave modes and quantify the fraction of energy around the MHD modes in linear and nonlinear simulations. Finally, the variation of the correlation time of the external forcing is studied in detail along with the effect on the particle energization (and clustering) and the presence of wave modes. More specifically, particle energization tends to decrease when the fraction of linear energy increases, supporting the idea that energization by structures is the dominant mechanism for particle energization instead of resonance with wave modes as suggested by Fermi energization theory.


INTRODUCTION
One of the first and most influential explanations for the origin of cosmic ray radiation was proposed by Enrico Fermi in 1949 (Fermi 1949).In that work, he argued about particles that could interact and resonate with passing Alfvén waves to achieve high kinetic energy.This explanation was later refined into the Quasi Linear Theory (QLT), which gave clear conditions for this resonance (Stix 1992).In terms of particle energization, QLT provides an estimation for diffusion coefficients in momentum space.However, these calculations rely on the assumption of weak turbulence (Chandran 2008(Chandran , 2005)), where each field is decomposed as a mean background value plus a small amplitude fluctuation described as a collection of weakly interacting waves.
For plasmas in a fully developed turbulent regime, as those present in astrophysical plasmas, these assumptions are not fulfilled and QLT is no longer adequate.Moreover, in the so-called strong turbulence regime, the main role in particle energization is played by self consistent structures (Dmitruk et al. 2004;Dmitruk & Matthaeus 2006;González et al. 2017;Lemoine 2021;Pugliese & Dmitruk 2022;Pezzi et al. 2022;Balzarini et al. 2022).Particles are able to exploit the electric fields present in these structures to obtain a net energization during a whole gyroperiod (Dmitruk et al. 2004;Pugliese & Dmitruk 2022).In the presence of a guide field, this energization is mainly perpendicular for protons and other heavy ions.Recent theoretical (Lemoine 2021) and numerical studies (Greco et al. 2014;Pezzi et al. 2022;Pugliese & Dmitruk 2022) have shown that some structures are also able to capture protons, which greatly enhances their energization capabilities.
At larger spatial scales, the magnetohydrodynamic (MHD) model provides a comprehensive explanation of interplanetary space plasmas, including the solar wind and planetary magnetospheres (Goedbloed et al. 2010).Within this framework, two distinct descriptions of the turbulence phenomenon can be considered.On one hand, the wave behavior, which describes field fluctuations as a collection of waves allowed by the MHD approximation, such as Alfvén waves or magneto-acoustic waves in the plain compressible case (see, Galtier 2006).On the other hand, the well-known strong turbulence regime, in which a broadband spectrum of fluctuations with coherent structures and the presence of intermittency, including non-propagating fluctuations, is present in wavenumber and frequency spaces (see, Pouquet & Yokoi 2022).To analyze the wave versus strong turbulence regimes, a spatiotemporal spectrum analysis of the fields can be performed (see, Meyrand et al. 2015;Clark di Leoni et al. 2015a).The scenario emerging from this analysis is one in which both the wave and turbulent behaviors of the fluctuations could coexist without contradictions (Dmitruk & Matthaeus 2009).This is supported by in situ observations data, such as those of the solar wind, which show both wave and turbulent behaviors in fluctuations (Barnes et al. 1979;C. Y. Tu 1995;Sahraoui et al. 2020).It is worth mentioning that the spatiotemporal analysis have been used to investigate incompressible hydrodynamic and incompressible and compressible MHD turbulence from experimental and numerical data (Clark di Leoni et al. 2014;Andrés et al. 2017;Lugones et al. 2016Lugones et al. , 2019;;Brodiano et al. 2021;Kawachi et al. 2022).
Charged test particles provide a useful representation of a small fraction of charged particles in a plasma that can be energized by the fields; however, neglecting the feedback of the particles into the electromagnetic field -which would require a kinetic description of the plasma -limits their application.Despite this, charged test particles are commonly used to study the energization of charged particles in different scenarios.In any case, electromagnetic fields need to be prescribed and it is usual to do so in Fourier space (Dalena et al. 2012;Tautz & Dosch 2013).Multiple models are available to provide the amplitudes simulating turbulent spectra such as Slab or 2D (Shalchi 2009), while phases are usually taken as random and uncorrelated.For example, in solar wind a combination of 20% Slab / 80% 2D is realistic at 1AU heliocentric distance (Bieber et al. 1996).Another approach is to use electromagnetic fields obtained from the evolution of MHD equations (Dmitruk et al. 2004;Dalena et al. 2014;González et al. 2016;Pugliese & Dmitruk 2022).
One important question is whether the wave or strong turbulent picture of the background plasma is more pertinent for the energization of the charged particles.In the present work, we concentrate specifically on this issue, analyzing the behavior of charged test particles with different techniques to determine the effect of the full MHD evolution of the field versus an evolution where non-linearities are artificially suppressed, resulting in only linear (wave) behavior.Our study provides important insights into the behavior of charged particles in plasmas and how they are energized by the surrounding fields.
The paper is organized as follows: in section 2, we introduce the theoretical CMHD set of equations and we present the linearized and dimensionless version.Then, we show the dispersion relation of the MHD waves modes, within a briefly explanation of the spatio-temporal spectrum technique.Finally, we present the equations regarding charged test particles and describe the Voronoi Tessellation method, used to quantify the concentration of particles in space.In section 3, we describe the numerical set up for non-linear and linear simulations along with some details about phase randomization for the linear runs and particle integration.In section 4, we present our results.Finally, in 5 we summarize our main findings.

Compressible MHD equations
The three dimensional (3D) compressible MHD (CMHD) model is given by the mass continuity equation, the induction equation for the magnetic field, the momentum Navier-Stokes equation and a polytropic state equation (relating pressure and densisty) and involves fluctuations of the velocity field u, density ρ and magnetic field B = B 0 + b, where B 0 is the mean field and b is the fluctuating part, Note that changes of the magnetic field absolute value B 0 involves modifying the relative amplitude between the initial fluctuations b and B 0 (with the initial conditions untouched).In addition, P is the scalar isotropic pressure, J = (4π/c)∇ × B is the electric current, γ = 5/3 is the polytropic index, and µ and η are the dynamic viscosity and magnetic diffusivity, respectively.The main purpose of these last terms is to dissipate energy at scales smaller than MHD scales, while allowing us to study with an adequate scale separation compressible effects at the largest scales.The set of equations ( 1)-( 4) can be written in a dimensionless form in terms of a characteristic length scale L 0 , a mean scalar density ρ 0 and pressure P 0 , and a typical magnetic and velocity field magnitude b rms and v 0 = b rms / √ 4πρ 0 (i.e., the r.m.s.Alfvén velocity), respectively.Then, the unit time is t 0 = L 0 /u rms , which for MHD becomes the Alfvén crossing time.The resulting dimensionless equations are, where M s = v 0 /C s is the sonic Mach number and C 2 s = γP 0 /ρ 0 is the sound speed.The kinetic and magnetic nominal Reynolds numbers are also defined as R e = L 0 v 0 /ν and R m = L 0 v 0 /η, respectively, with ν = µ/ρ 0 the kinematic viscosity.Considering a static equilibrium (u 0 = 0) with homogeneous external magnetic field B 0 = B 0 ẑ, a constant density ρ 0 , and a constant pressure P 0 , we can linearize Eqs. ( 5)-( 8), where the polytropic equation ( 4) was used to replace the pressure in the equations.

Compressible MHD waves modes
In order to study the CMHD normal modes, we used Eqs.( 9)-( 11) and obtain the dispersion relation ω(k) of small amplitude waves propagation in a plasma.It is straightforward to show that there are three independent propagating modes (or waves), which correspond to the so-called Alfvén waves (A), fast (F), and slow (S) magnetosonic waves (e.g., Fitzpatrick 2014), where β = (C s /u A ) 2 is the plasma beta, i.e., the ratio of plasma pressure to magnetic pressure, which can expressed as β = 1/(M s B 0 ) 2 ) with u A = B 0 / √ 4πρ 0 the Alfvén velocity, C s as defined above and k = |k| = k 2 + k 2 ⊥ with and ⊥ the wavenumber component along and perpendicular to the external magnetic field, respectively.It is worth noting that B 0 is a parameter and has no relation with the initial fluctuations of the system.On one hand, Alfvén waves are incompressible fluctuations transverse to the external magnetic guide field B 0 .On the other hand, both fast and slow modes, unlike Alfvén modes, carry density fluctuations and their magnetic field perturbations have longitudinal and transverse components.In the case of fast modes, the magnetic field and the plasma are compressed by the wave motion, such that the restoring force is large and hence the frequency and the propagation speed are high.While, for slow modes, the magnetic field oscillation and the pressure variation are anti-correlated with each other such that the restoring force acting on the medium is weaker than that for the fast mode.For this reason the frequency and the propagation speed are the lowest among the three MHD waves branches.Note that for the perpendicular propagation (i.e., k = 0 and k ⊥ = 0), the Alfvén and slow modes become non-propagating modes (i.e., ω A,S = 0) and are degenerate, but they can be distinguished using their different polarization, since δB A = 0 and δB S = 0 (Andrés et al. 2017).Finally, it is worth mentioning that, we adopt the assumption that energy concentrated closely to the linear dispersion relation can be explained by linear and weak turbulence theories, while any spread away from those linear curves is a sign of strong turbulence that requires fully nonlinear theories to be understood (see, Andrés et al. 2017;Brodiano et al. 2021).

Spatio-temporal spectrum
The spatio-temporal spectrum consists of calculating the complete spectrum in wavenumber and frequency for all available Fourier modes in a numerical simulation or an experiment (Clark di Leoni et al. 2014Leoni et al. , 2015a)).As a result, it can distinguish between modes that satisfy a given dispersion relation (and are thus associated with waves) from those associated with nonlinear structures or turbulent eddies, and quantify the amount of energy carried by each of them at different spatial and temporal scales.It is worth mentioning, that this method does not require the pre-existence of wave modes or eddies.In the following, the spatio-temporal magnetic energy spectral density tensor is defined as: where Bi (k, ω) is the Fourier transform in space and time of the i-component of the magnetic field B(x, t) and the asterisk implies the complex conjugate.The magnetic energy is associated with the trace of E ij (k, ω).
As the external magnetic field B 0 in the simulations points in ẑ, in practice, we will consider either i = j = y or i = j = z, to identify different waves based on their polarization (either transverse or longitudinal with respect to the guide field).In all cases, the acquisition frequency was at least two times larger than the frequency of the fastest wave, and the total time of acquisition was larger than the period of the slowest wave in the system.It is worth mentioning that spatio-temporal spectra have been used before in numerical simulations and experiments of rotating turbulence (Clark di Leoni et al. 2014), stratified turbulence (Clark di Leoni & Mininni 2015), quantum turbulence (Clark di Leoni et al. 2015b), and IMHD turbulence simulations (Meyrand et al. 2015(Meyrand et al. , 2016;;Lugones et al. 2016), compressible MHD turbulence (Andrés et al. 2017;Brodiano et al. 2021) and in spacecraft observations (Sahraoui et al. 2003(Sahraoui et al. , 2010)).Quantifying the relative presence of waves and/or nonlinear structures is the main outcome expected from the spatio-temporal analysis.In particular, we use the spatio-temporal spectra to search for the presence of waves and quantify the energy found near their respective dispersion relation in a linear and nonlinear MHD turbulence.

Test particle equations
We studied charged particles evolving with the dynamical MHD fields, but those fields are unaffected by the particles (i.e., test particles).The dynamics of a test particle in an electromagnetic field are given by the non-relativistic equation of motion: where the electric field E is obtained from Ohm's Law and its dimensionless (using a characteristic electric field While a more general version of this law includes the electronic pressure (i.e., ∇P e /ρ) and the Hall term (i.e, J × B/ρ), here we neglect them in order to simplify the analysis and interpretation of the results and to maintain consistency with the compressible MHD equations introduced in the previous sections.Several works have studied the effect of these terms on test particles (Dmitruk & Matthaeus 2006;Pugliese & Dmitruk 2022;Balzarini et al. 2022).As such, density fluctuations will not directly affect particle motion through Eq. ( 16), but could allow the existence of magnetosonic wave modes for the particles to interact with.
The parameter α is related to the charge-to-mass ratio and represents the gyrofrequency ω g in a magnetic field of intensity b rms (Pugliese & Dmitruk 2022).Its inverse value 1/α represents the nominal gyroradius r g (in units of L 0 ) for particles with velocity v 0 in a magnetic field b rms , and gives a measurement of the range of scales involved in the system (from the outer scale of turbulence to the particle gyroradius).For plasmas with a strong magnetic field amplitude B 0 , we can define an average gyroperiod τ p = 2π/ω g = 2π/αB 0 .In the case of protons, we can relate α to the MHD field parameter through, where d i = m p c/ 4πρ 0 e 2 is the proton inertial length, with m and e the proton mass and charge, respectively.Furthermore, in the present paper, we identify the proton inertial length with the dissipation scale d i = l d , given the solar wind observations supporting d i ∼ l d (e.g., Leamon et al. 1998).

Particle energization and clustering
Recently, Pugliese & Dmitruk (2022) reported a link between high energization of test particles and a preferential concentration (clustering) for test particles.In the present work, to quantify this preferential concentration, we employ the Voronoi tessellation method and compare volume statistics against a Random Poisson Process (RPP) (see, Monchaux et al. 2010;Obligado et al. 2014;Reartes & Mininni 2021).In particular, this process provides a cell for each particle, whose volume V i can be interpreted as the inverse of the local particle density.Therefore, by studying the statistics of these volumes, we can quantify any preferential concentration that may be present in the system.In addition, we compared against those of a uniform distribution in space, known as RPP (Uhlmann 2020).One of these statistic tools is the standard deviation of the volumes σ V , which should coincide with σ RP P for a uniform distribution and increases as the preferential concentration increases.
In the aforementioned work it was found that particles accumulated in regions where 16) this also implies (∇ × E) z < 0 and thus a clockwise rotating electric field.Therefore, in this regions the electric force and particle gyration are aligned and there is positive energization after a whole gyroperiod.This, in combination with the trapping effect of ∇ ⊥ • u ⊥ < 0, makes this mechanism very effective at energizing particles.

Nonlinear simulations
The nonlinear simulations are performed by numerically solving Eqs. ( 1)-( 4) using a pseudospectral method with periodic boundary conditions in a cube of size L box = 2π.This scheme ensures exact energy conservation for the continuous time spatially discrete equations (Mininni et al. 2011), where we use a spatial resolution of N 3 = 512 3 Fourier modes.Time integration is achieved through a second order Runge-Kutta method.Aliasing is removed using the two-thirds truncation method (Orszag 1971), such that the maximum wavenumber resolved is κ ≡ N/3 = 170.To ensure the resolution of the smallest scales (κ > k d ), we use kinematic and magnetic Reynolds numbers of R e = R m ≈ 2400.
In order to reach a steady turbulent state, the system is forced using mechanical and electromotive forces f and ∇ × m, respectively.These forces are generated with random phases in the Fourier k-shells 2 ≤ |k| ≤ 3 every correlation time τ c , which is also a controlled parameter in the simulations.Forces at intermediate times are obtained by linearly interpolating the previous and next steps.We repeat this procedure for three different values of τ c , yielding the three different stationary states summarized in Table 1.The length and time scales are defined individually for each simulation.For the characteristic length scale L 0 , we use the energy containing scale L 0 = 2π (E(k)/k)dk/ E(k)dk where E(k) is the isotropic energy spectrum.For the velocity scale, we use the Alfvén velocity of the magnetic field fluctuations v 0 = b rms / √ 4πρ 0 .For the time scales, we alternatively use the large eddy turnover time t 0 = L 0 /v 0 and the particle gyroperiod τ p = 2π/αB 0 , depending on whether we analyze field or particle properties.Finally, as seen in Table 1, we have a ratio of mean magnetic field to magnetic field fluctuation (B 0 /b rms ) of 1 : 9, similar to values reported in the solar wind (Hadid et al. 2017;Andrés et al. 2022).

Linear simulations
The linear simulations are performed in the same way as the nonlinear ones but we cancel all the non-linear terms in the MHD equations, as shown in Eqs. ( 9)-( 11).In the absence of nonlinear terms, there is no energy cascade from the injection scale L 0 to smaller scales.Therefore, modes with wavenumber k outside the forced shell would quickly vanish in the presence of the dissipation terms in Eqs. ( 10)-(11).To counter this, we evolve the fields in Eqs. ( 9)-( 11) without forcing and without dissipation (i.e.R e , R m → ∞), thus ensuring a constant energy spectra.Furthermore, we use a fourth order Runge-Kutta scheme for temporal integration, as the second order scheme becomes unstable in the absence of dissipation terms.
The initial conditions for the linear simulations are given by different variations of the initial conditions used in NL1, as summarized in Fig. 1.For the L run, we use these initial conditions unchanged, while for the LR run we perform a phase randomization.We achieve this by transforming each Fourier coefficient ψ k → e iφ k ψ k for all the scalar fields ρ, u j and b j , where φ k are random phases independently chosen for each k and ψ.We then enforce the necessary conditions on the resulting fields (i.e.hermiticity, Coulomb gauge).These first two runs are built to have the same energy spectra as NL1, but thoroughly destroy any cross-field correlation and structures present in the initial conditions of NL1 for the LR run (Alexakis et al. 2007).
The final two runs LR80 and LR40 use the same Fourier coefficients as LR, but imposing ψ k = 0 for |k| > 80 and |k| > 40, respectively.As this process slightly reduces the total energy, we compensate by uniformly re-scaling all the Fourier coefficients, thus preserving any power-law behaviour.

Particle integration
Initially, particles were uniformly distributed in the box and with a Gaussian velocity distribution function, with a root mean square (rms) of v 2 i 1/2 ≈ 1.2v 0 .As shown in Table 1, this value lies between the Alfvén velocity v 0 and the rms value of the velocity field u rms .Furthermore, we chose α = 60 for the particles, which is consistent with the values of L 0 /l d for all simulations.The numerical integration of Eq. ( 15) was done by a Runge-Kutta method with the same order as that used for field integration (see above).The values of the fields at each particle position are obtained by cubic splines in space from the three-dimensional grid of the MHD simulation.The particle trajectories were integrated for ∼ 300τ p .

Linear vs nonlinear
In this section, we compare the dynamics of test particles in the nonlinear run NL1 versus the dynamics in the linear counterparts L, LR, LR80 and LR40.We begin with the evolution of the mean particle energy, showing separately the parallel and perpendicular component to the magnetic guide field in Fig. 2(a) and Fig. 2(b), respectively.As it is known for protons, energization is much higher in the perpendicular direction (represented here by the x component) with respect to the parallel direction (represented here by the z direction).For both directions, it is clear that the linear case L has significantly smaller energization when compared to the full nonlinear simulation NL1.However, for the phase randomizing (LR run), the energization increases greatly and surpasses that of NL1 run.Truncating the spectra at |k| = 80 (LR80 run) has very little effect, but truncating at |k| = 40 (LR40 run) greatly reduces the energization, obtaining values similar to the linear L run.Furthermore, in the perpendicular case shown in Fig. 2(a) all the linear runs seem to display a diffusive or subdiffusive behaviour (in momentum space) v 2 x ∼ t α with α 1.This fact suggests that the underlying mechanism for energization (in the linear case) is analogous to that of Brownian motion, with delta correlated energy increments.In this analogy, particles would have a very strong and time localized interaction with the fields.In the context of QLT, this very strong interaction could be understood as a resonance with some specific wave present in the system.This high energization quickly removes the particle from the resonance condition within a few gyroperiods, thus ensuring the interaction to be localized in time.
On the other hand, energization in the nonlinear case is superdiffusive, which is related to interaction between particles and self-consistent structures present in the plasma, as discussed in Section 2.5.With this in mind, we could relate the drop in energization as NL1→L to the disappearance of structures in the plasma due to the linear evolution.In Eqs. ( 12) and ( 13) we see that waves in the CMHD model are indeed dispersive and as such they will destroy any structure initially present in the system.To investigate this, we calculated the radial two-point autocorrelation Γ of ∇ ⊥ • u ⊥ in the plane perpendicular to the magnetic guide field.In Fig. 3(a) we show the autocorrelation at multiple times (darker colors represent later times).From these curves, we calculated the correlation length c , defined in this case as the displacement at which the autocorrelation drops below 10% (horizontal dotted black line).In Fig. 3(b) we show the autocorrelation length c as a function of time, where we note that the NL1 run holds c ≈ 18l d during the whole simulation.On the other hand, the linear run L starts with a similar value to the NL1 run but very quickly drops to c ≈ 8l d and then slowly decreases toward that of the LR run ( c ≈ 2l d ), which by construction should have negligible correlation (see Sec. 3.2).Therefore, this confirms that the linear evolution eliminates most structures initially present in the field in less than 20τ p , preventing particles from exploiting the energization mechanism described in Sec.2.5.To further check this claim, in Fig. 4 we show the standard deviations of the Voronoi volumes as a function of time.We see that the NL1 run quickly increases its clustering while all the linear runs remain very close to uniformity.The slight increase in their σ V could mean that this mechanism still survives in some small scale regions, but given the diffusive energization of Fig. 2 we could consider this effect secondary at best.We turn now to the comparison of the linear energization runs between themselves.Having discarded structure interaction in the linear cases, we are left only with wave-particle resonant interaction.The first observation is that the phase randomization L→LR seems to vastly increase particle energization.This fact mainly shows the importance of the phase mixing hypothesis in QLT.The L run clearly does not fulfill this hypothesis, as its initial conditions are derived from a fully nonlinear turbulent simulation and as such display high phase correlation.Alternatively, we could visualize the wave-particle interaction as a conjunction of resonance in terms of frequency and an initial alignment between the field and particle velocity.This alignment is more easily achieved under the phase mixing hypothesis, while for a turbulent state phases do not actually fill every possible value at all frequencies.
We can visualize these different energization mechanisms by computing the two-dimensional PDF of particle velocities at the end of each simulation.The clear geometrical difference of each PDF shown in Fig. 5 suggests that each mechanism is fundamentally different, as the data is normalized to the unit variance.Close to the origin, all distributions are elliptic (although with different foci).As we move from the origin, each distribution changes significantly.In particular, for the NL1 run, they tend to rhombi reminiscent of constant 1-norm (i.e., |v x | + |v z | =const.),suggesting that simultaneously high perpendicular and parallel energization is unlikely, as expected from this structure energization (high v z particles are difficult to trap).For the L run, the shapes are clearly circular, corresponding to constant 2-norm (v 2 x + v 2 z =const.)or kinetic energy.This is consistent with two independent Gaussian variables, as expected from a diffusive process in momentum space.Finally, the LR run is less clear, with shapes similar to rectangles reminiscent of constant ∞-norm (max{v x , v y } =const.).
By comparing the LR run with its truncated counterparts LR80 and LR40, we first see that the first truncation (i.e., LR→LR80) has very little effect on particle energy, implying that particles resonate mainly with waves of |k| < 80. Furthermore, the second truncation LR80→LR40 reduces energization to even lower than the L run, showing that particles resonate mainly with waves of |k| > 40.Two-dimensional velocity distributions for the LR80 and LR40 runs (not shown here) are very similar to those of LR and L in Fig. 5, respectively.While the similarity between LR80 and LR is to be expected, that between LR40 and L is not obvious and seems to show that both runs share the same energization mechanism.This last fact implies that changes in the phase distribution of the waves could be as important as changes in the energy spectrum.
In order to distinguish between the presence of waves and structures in a plasma and to investigate which dominate the dynamics, the spatio-temporal spectra were used (e.g., Clark di Leoni et al. 2015a; Figure 6.Spatio-temporal spectrum E zz (k x = 0, k y , k z = 0) for the magnetic field fluctuations parallel to B 0 , for the runs NL1 (left) and L (right).The spectrum is shown as a function of ω and k y for fixed k x = k = 0.The green solid and the green dashed-dotted lines correspond to the linear dispersion relation of fast magnetosonic waves (ω F ) and of slow magnetosonic waves (ω S ), respectively.We include the particle gyrofrequency ω g (dashed line) and the sweeping frequency w sw (blue line).Andrés et al. 2017).Therefore, we resolved the waves in time and space by choosing a very high frequency cadence to store the magnetic field, in particular we used dt = 2.5 × 10 −4 as the temporal sampling rate.It is worth to mention that, we assumed that the energy concentrated around the linear dispersion relation can be explained by linear and weak turbulence theories (Chandran 2005(Chandran , 2008)), while any spread away from the dispersion relation is a sign of strong turbulence that requires fully nonlinear theories to be understood.Figure 6 shows the spatio-temporal spectrum of the parallel magnetic field E zz (k x = 0, k y , k z = 0, ω) for the NL1 and L runs (the rest of linear runs have similar spectra to the L run), i.e. we separate the linear from the nonlinear case.The dispersion relation for the fast and slow magnetosonic waves given by Eq. ( 13) is shown in green solid and green dashed-dotted lines, respectively.Also, we include the particle gyrofrequency ω g which lies between wavenumbers 40 < |k| < 80 and the characteristic sweeping frequency given by ω sw ∼ U rms k 2 ⊥ + k 2 (see, Lugones et al. 2016).As we expected, for the linear run L, the magnetic energy is located around the magnetosonic branches, while for the nonlinear run NL1, the energy accumulates mainly in the slow waves (with ω/B 0 = 0 or non-propagating modes) for all wavenumbers.We also noted a small portion of energy spread along the fast branch for low wavenumbers.However, in the nonlinear case, most of the energy is spread across the spectrum due to its turbulent dynamic.
Furthermore, the absence of sweeping in the linear case is also useful to understand the low energization in the L run.As shown in Fig. 3, correlation quickly drops during linear evolution but not quite enough to reach the LR case, which implies that some structures may survive longer.The absence of sweeping produced by the linear evolution means these surviving structures are not advected by the flow.However, test particles are advected, thus making trapping more difficult and preventing any surviving structure from effectively energizing particles.

The effect of τ c on particle energization
For this section, we repeated the analysis from the last section with the simulations NL2 and NL3, which differ from the NL1 run in the correlation time of the forcing τ c (see Table 1).Energization is qualitatively similar to that of NL1 (as shown in Fig. 2), but not quantitatively.The same occurs for the correlation length and volume deviation σ V , as summarized in Table 2 (energization and σ V are calculated at the end of the simulation).We see that both energization and correlation length decrease along with τ c , while σ V displays no clear tendency and remains mostly similar.First, the reduction of c /l d could be expected, as faster forcing (lower correlation time) prevent structures from being stable, either for lack of size or intensity.Secondly, we could conclude that particles tend to accumulate similarly in all cases, but their ability to exploit energization mechanism diminishes for lower τ c .
Therefore, we need another way to quantify particle interaction with structures, including some dynamical information.For this purpose, we performed Voronoi tessellation for each time step and calculated whether the particle was clustered.We define a particle as clustered when its cell volume V is lower than some threshold value V th , which is obtained by comparing the cell volume PDF with that of a RPP (see references in Section 2.5 for more details).As seen in Fig. 4, clustering takes some time to settle, and as such the label clustered may not mean much initially.Towards the end of the simulation, clustered is almost equivalent to trapped inside a structure.Then, we can determine when and where are particles accumulating and therefore calculate the amount of consecutive time they spent clustered.We define a streak as an interval of (discrete) time for which the particle is clustered and its duration as S. In particular, we calculated the mean streak time S for each simulation by averaging over all streaks from all particles and display it in the last column of Table 2. Particles are clustered for around 2 full gyroperiods in all cases, but the mean time decreases with τ c .This shows that while the clustering is very similar instantaneously for all runs, there is more exchange (between clustered and not clustered particles) in low τ c simulations, suggesting clusters become more feeble and unstable with lower correlation time.For lower values of τ c , the rapidly changing forcing may be experienced by particles as kicks that remove them from the structures trapping them.
To further perform this streak analysis, we related it directly with particle energization.For this, we separated our data into intervals of duration time nτ p (n = 1, 2, 3) and we selected particles that are clustered during this time (whose streaks contained the interval under study).For these particles, we calculated their mean perpendicular energization during these intervals and display it in Fig. 7.We can see that after ∼ 100τ p (long enough for clustering to settle), energization becomes exponential in time.
In order to understand this, we propose a very simple model for the interaction with these structures.If we write the particle velocity in terms of its parallel and perpendicular components v = v ⊥ + v , we can define the perpendicular energy as ε ⊥ = |v ⊥ | 2 /2.Using Eq. ( 15), we could derive an evolution equation for the perpendicular energy ε ⊥ , The first term P ⊥ = αE • v ⊥ corresponds to net energization while the second relates to exchange with parallel energy, such as pitch angle scattering.Averaging over one gyroperiod, where we have approximated the trajectory on the perpendicular plane as circular and then used Stokes' Theorem, with a minus sign accounting for the clockwise circulation.For a strong guide field, we could approximate Eq. ( 16) as E ⊥ ≈ −u ⊥ × B 0 and therefore the parallel component of its curl is Expanding the integral at the lowest order in the gyroradius R g we achieve, Using that R g = |v ⊥ |/αB 0 and recalling the definition of ε ⊥ , we can substitute in Eq. ( 18) to obtain, which for approximately constant ∇ ⊥ • u ⊥ < 0 predicts exponential increase for ε ⊥ .
The slopes in Fig. 7 represent λ and are very similar, suggesting that clustered particles are energized at the same rate.The difference in energization observed in Table 2 must therefore be related to the time each particle spends clustered, as shown by S .Particles can leave structures in basically 3 ways: (a) escaping vertically due to high v z , (b) reaching the maximum allowed gyroradius, or (c) being pushed out by some fluctuation.We disregard option (a) as particles in all simulations have very similar parallel energization (not shown).Option (b) implies that the gyroradius R g = v ⊥ /αB 0 is comparable to the size of the structures, which could be taken as c and would be reasonable considering how it depends on τ c (see Table 2).We can calculate the required kinetic energy as v 2 ⊥ ∼ (αB 0 c ) 2 ∼ 1600v 2 0 , which is achieved by practically no particle (less than 1 particle in 10 4 ) and as such cannot be the dominant cause.This leaves option (c) as the last and main possibility, suggesting that structures are less robust against fluctuations and thus weaker.In the wave/turbulence dichotomy, we could attribute this weakness to the prevalence of waves over structures in the system, to which we will dedicate the next section.

The effect of τ c on spatio-temporal spectra
For the study of the relevance of waves in the system for different τ c , we quantitatively analyze the spatio-temporal magnetic spectra.Figure 8 shows the spatio-temporal spectrum of the perpendicular magnetic field fluctuation component E xx (k x = 0, k y = 15, k z , ω) for the NL1, NL2 and NL3 runs.For an easy comparison, we include the dispersion relation for the Alfvén, fast and slow magnetosonic waves in blue dashed, green dashed-dotted and orange solid lines, respectively.We also added the particle gyrofrequency and the sweeping frequency in red dashed and solid blue lines, respectively.We observed that the energy is mainly located around the slow branch and, in lesser extent, around Alfvén waves.Moreover, as τ c decreases (i.e., as we move from NL1 to NL3), the energy around wave modes increases for lower values of k z .In fact, it is worth noting that as we increase the correlation time the energy tends to spread slightly towards the sweeping frequency.This result is indeed compatible with the behaviour on c /l d , as the sweeping energy is mostly related to the non-linear structures (see Table 2).
In order to analyze the magnetosonic modes, we studied the spatio-temporal spectrum of parallel magnetic field fluctuation E zz (k x = 0, k y , k z = 15).Figure 9 shows the same trend that we observed in Fig. 8.In fact, the energy around the fast branch decreases noticeably as we increase τ c .It is worth mentioning that the slow branch is completely immerse in the sweeping region, therefore, we can not come to a conclusion about the increase (or decrease) of the magnetic power energy around the slow mode.However, in the case of NL1, the energy is more concentrated around ω = 0 compared to lower values of the correlation time and the sweeping frequency correctly delimits this energy.For a more precise study, we used an integration method to quantify the amount of energy near the different wave modes in each spatio-temporal spectrum.Clark di Leoni et al. (2015a) calculated the ratio of the energy accumulated near these wave modes to the total energy in the same wavenumber as, with ω A,F,S the frequencies that satisfy a certain dispersion relation (and the E xx (k z ) component is an illustration example only).Figure 10 shows the energy around (a) Alfvén and (b) slow magnetosonic waves for runs NL1, NL2 and NL3.In particular, we observed that the amount of energy near the slow branch is similar as we increase the correlation time.However, for small scales (k z ∼ 100), the energy in run NL1 tends to be smaller than the rest of the simulations as a result of the energy moving towards the sweeping frequency.In the case of Alfvén modes, most of the energy is concentrated around low wavenumbers and tends to decrease for higher k z .It is worth noting that we are not including the fast branch since most of the energy is not located around this wave mode.Figure 11 shows (a) the energy near fast magnetosonic waves for runs NL1,NL2 and NL3 and (b) the difference between NL3 and the first two non-linear simulations.We observe that the energy located around the fast branch reaches its maximum value for the smallest wavenumbers.Moreover, we obtained a significant growth of the energy as the correlation time decreases.Finally, the NL1 run differs from the NL3 run by 11%.

CONCLUSIONS
In the present work, we investigated the interplay of linear waves and coherent structures in compressible MHD turbulence and how this affects test particle (proton) energization.We compared a nonlinear evolution against a linear evolution of the same initial turbulent state.This initial state is obtained from a fully nonlinear evolution, which is more realistic than the usual spectra models with random phases.We found that under this initial state, energization in the linear case is much lower than the nonlinear case.However, this situation inverts after the phases of the initial state are randomized, showing the relevance of phase distribution in particle energization.This last numerical result allow us to reinterpret the role of phases in particle energization, which are usually treated as secondary with respect to the spectra.Additionally, this result shows that linear evolution, no matter how realistic the spectra is, can not faithfully reproduce energization observed from a nonlinear evolution.
We showed that wave-particle interactions are mainly in modes 40 < |k| < 80, which for Alfvén waves is consistent with ω A ≈ ω g .This could be related to second order Fermi energization ω±k z v z = ω g in the reasonable limit ω g k z v z .First order Fermi energization requires |v z | ≈ v A ≈ 9v 0 , which we observed from the parallel energization in Fig. 2 is not likely.A similar reasoning follows for the fast magnetosonic waves, whose frequency ω and resonant velocity ω/k z are even higher.On the other hand, first order Fermi energization is more likely to happen in the slow magnetosonic waves and is the only possible mechanism left for the LR40 run.
The fact that the L and LR40 runs have similar energization mechanisms and rates is striking, as it further reinforces the role of phases in particle energization.If slow magnetosonic waves are responsible for the energization in the LR40 run, they must also be for that of the L run.However, this is unlikely, as their phase and energy distributions are different in each run because phase randomization decorrelates and tends to evenly distribute energy between the different branches.According to our numerical results, it would seem that phase distribution is very important for high wavelength waves and becomes practically irrelevant for low wavelength waves.This result suggests that the simple picture of particles resonating with one wave at a time may be misguiding and perhaps resonance broadening is enhanced for such complex spectra.The PDFs of Fig. 5 seem to imply that the energization of L/LR40 and LR/LR80 are fundamentally different.
Furthermore, we analysed the effect of the correlation time of the external forcings τ c on the spatio-temporal spectra and the particle energization.Here, we further discussed the trapping and energization mechanism of the structures, showing that clustered particles energize exponentially.

Figure 1 .
Figure 1.Schematic diagram showing how the initial conditions of each linear run are constructed from the initial conditions of the NL1 run.

Figure 2 .
Figure2.Time evolution of the mean kinetic energy for the perpendicular (left) and parallel (right) components, for the non-linear and linear simulations.For the time scale, we used the particle gyroperiod τ p .

Figure 3 .
Figure 3. (a) Radial two-point autocorrelation Γ of ∇ ⊥ • u ⊥ in the plane perpendicular to the magnetic guide field at multiple times (later times are represented by darker colors) as a function of the displacement length normalized by the dissipation scale d , for non-linear and linear simulations.(b) Correlation length c , normalized by d , as a function of time; the correlation length c is obtained as the displacement length at which the autocorrelation Γ drops below 10% of its maximum value (indicated by the dotted line in part (a) of this figure) .

Figure 4 .
Figure 4. Standard deviation of the Voronoi volumes normalized by the standard deviation for a uniform distribution, as a function of time for non-linear and linear simulations.

Figure 5 .
Figure 5. Two-dimensional PDF of particle velocities normalized by its standard deviation at the end of the simulation for the NL1, L and LR runs, respectively.The contours display different geometry, implying different underlying mechanism for energization.

Figure 7 .
Figure7.Mean of the structure function of the perpendicular particle energization for intervals with time duration nτ p (where n = 1, 2, 3).We selected particles that are clustered during this time.

Figure 8 .
Figure8.Spatio-temporal spectrum E xx (k x = 0, k y = 15, k z ) for the magnetic field fluctuations perpendicular to B 0 .The spectrum is shown as a function of ω and k (k z ) for fixed k x = 0 and ky = 15.The blue dashed, green dashed-dotted and orange solid lines correspond to the linear dispersion relation of Alfvén waves (ω A ), of slow magnetosonic waves (ω F ) and of fast magnetosonic waves (ω S ), respectively.We include the particle gyrofrequency ω g and the sweeping frequency in red dashed and blue solid lines, respectively.

Figure 9 .
Figure9.Spatio-temporal spectrum E zz (k x = 0, k y , k z = 15) for the magnetic field fluctuations parallel to B 0 , for the runs L and NL1.The spectrum is shown as a function of ω and k y for fixed k x = 0 and k = 15.The green solid and the green dashed-dotted lines correspond to the linear dispersion relation of fast magnetosonic waves (ω F ) and of slow magnetosonic waves (ω S ), respectively.We include the particle gyrofrequency ω g and the sweeping frequency in red dashed and blue solid lines, respectively.

Figure 10 .
Figure 10.Quantification of the amount of energy around (a) Alfvén and (b) slow magnetosonic branches presented in spatio-temporal spectrum E xx (k x = 0, k y = 15, k z ).

Figure 11 .
Figure 11.Quantification of the amount of energy around (a) fast magnetosonic branch presented in spatio-temporal spectrum E zz (k x = 0, k y , k z = 15).Difference between the NL3 run and the rest of the nonlinear simulations (b).

Table 1 .
Global quantities for nonlinear simulations Run τ c /τ p L box /L 0 L 0 /l d B 0 /b rms u rms /v 0

Table 2 .
Particle related quantities in nonlinear simulations