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

The three-dimensional 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 ₀ + b , where B ₀ is the mean field and b is the fluctuating part:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times ({\boldsymbol{u}}\times {\boldsymbol{B}})+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{u}} & = & -\displaystyle \frac{{\boldsymbol{\nabla }}P}{\rho }+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{4\pi \rho }\\ & + & \displaystyle \frac{\mu }{\rho }\left[{{\rm{\nabla }}}^{2}{\boldsymbol{u}}+\displaystyle \frac{{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}})}{3}\right],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{P}{{\rho }^{\gamma }}=\mathrm{constant}.\end{eqnarray} \tag{ 4 }$$

Note that changes in the absolute value of the magnetic field B₀ involve modifying the relative amplitude between the initial fluctuations b and B₀ (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 compressible effects at the largest scales with an adequate scale separation.

The set of Equations (1)–(4) can be written in a dimensionless form in terms of a characteristic length scale L₀, a mean scalar density ρ₀ and pressure P₀, and a typical magnitude of magnetic and velocity fields b_rms and ${v}_{0}={b}_{\mathrm{rms}}/\sqrt{4\pi {\rho }_{0}}$ (i.e., the rms Alfvén velocity), respectively. Then, the unit time is t₀ = L₀/u_rms, which for MHD becomes the Alfvén crossing time. The resulting dimensionless equations are

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{u}}=-\displaystyle \frac{1}{\gamma {M}_{s}^{2}}\displaystyle \frac{{\boldsymbol{\nabla }}P}{\rho }+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{\rho }\\ \qquad +\ \displaystyle \frac{1}{\rho {R}_{e}}\left[{{\rm{\nabla }}}^{2}{\boldsymbol{u}}+\displaystyle \frac{{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}})}{3}\right],\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times ({\boldsymbol{u}}\times {\boldsymbol{B}})+\displaystyle \frac{1}{{R}_{m}}{{\rm{\nabla }}}^{2}{\boldsymbol{B}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&P={\rho }^{\gamma }\end{eqnarray} \tag{ 8 }$$

where M_s = v₀/C_s is the sonic mach number and ${C}_{s}^{2}=\gamma {P}_{0}/{\rho }_{0}$ is the sound speed. The kinetic and magnetic nominal Reynolds numbers are also defined as R_e = L₀ v₀/ν and R_m = L₀ v₀/η, respectively, with ν = μ/ρ₀ the kinematic viscosity.

In this work, we are interested in the weakly compressible case, which we achieve by setting M_s = 0.25. This ensures small density fluctuations, which in our simulations implies $\langle {(\rho -{\rho }_{0})}^{2}{\rangle }^{1/2}/{\rho }_{0}$ between 13% and 17%, comparable to solar wind fluctuations (Hadid et al. 2017; Andrés et al. 2021; Brodiano et al. 2023).

Considering a static equilibrium (u₀ = 0) with homogeneous external magnetic field B ₀ = B₀ $\hat{z}$ a constant density ρ₀, and a constant pressure P₀, we can linearize Equations (5)–(8):

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}=-\displaystyle \frac{1}{{M}_{s}^{2}}{\boldsymbol{\nabla }}\rho +{\boldsymbol{J}}\times {{\boldsymbol{B}}}_{0}+\displaystyle \frac{1}{{\rho }_{0}{R}_{e}}\left[{{\rm{\nabla }}}^{2}{\boldsymbol{u}}+\displaystyle \frac{{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}})}{3}\right],\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial t}={\boldsymbol{\nabla }}\times ({\boldsymbol{u}}\times {{\boldsymbol{B}}}_{0})+\displaystyle \frac{1}{{R}_{m}}{{\rm{\nabla }}}^{2}{\boldsymbol{b}},\end{eqnarray} \tag{ 11 }$$

where the polytropic Equation (4) was used to replace the pressure in the equations.

In order to study the CMHD normal modes, we use Equations (9)–(11) and obtain the dispersion relation ω( k ) of small-amplitude waves propagating 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) and fast (F) and slow (S) magnetosonic waves (e.g., Fitzpatrick 2014):

$$\begin{eqnarray}&&{\omega }_{{\rm{A}}}^{2}(k)={k}_{\parallel }^{2}{u}_{{\rm{A}}}^{2},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\omega }_{{\rm{F}},{\rm{S}}}^{2}(k)={k}^{2}{u}_{{\rm{A}}}^{2}\left[\displaystyle \frac{1+\beta }{2}\pm \sqrt{\tfrac{{\left(1+\beta \right)}^{2}}{4}-\beta {\left(\displaystyle \frac{{k}_{\parallel }}{k}\right)}^{2}}\right],\end{eqnarray} \tag{ 13 }$$

where $\beta ={\left({C}_{s}/{u}_{{\rm{A}}}\right)}^{2}$ is the plasma beta, i.e., the ratio of plasma pressure to magnetic pressure, which can expressed as $\beta =1/{\left({M}_{s}{B}_{0}\right)}^{2}$ with ${u}_{{\rm{A}}}={B}_{0}/\sqrt{4\pi {\rho }_{0}}$ the Alfvén velocity, C_s as defined above, and $k=| {\boldsymbol{k}}| =\sqrt{{{\boldsymbol{k}}}_{\parallel }^{2}+{{\boldsymbol{k}}}_{\perp }^{2}}$ with ∥ and ⊥ the wavenumber components along and perpendicular to the external magnetic field, respectively. It is worth noting that B₀ is a parameter and has no relation to the initial fluctuations of the system. On one hand, Alfvén waves are incompressible fluctuations transverse to the external magnetic guide field B ₀. 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. For slow modes, the magnetic field oscillation and the pressure variation are anticorrelated 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 wave branches. Note that for the perpendicular propagation (i.e., k_∥ = 0 and k_⊥ ≠ 0), the Alfvén and slow modes become nonpropagating modes (i.e., ω_A,S = 0) and are degenerate, but they can be distinguished using their different polarizations, 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 close 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).

Determining the spatiotemporal 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. 2014, 2015a). As a result, it can distinguish 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 preexistence of wave modes or eddies. In the following, the spatiotemporal magnetic energy spectrum is defined as

$$\begin{eqnarray}&&{E}_{{ij}}({\boldsymbol{k}},\omega )=\displaystyle \frac{1}{2}{\hat{{\boldsymbol{B}}}}_{i}^{* }({\boldsymbol{k}},\omega ){\hat{{\boldsymbol{B}}}}_{j}({\boldsymbol{k}},\omega ),\end{eqnarray} \tag{ 14 }$$

where ${\hat{{\boldsymbol{B}}}}_{i}({\boldsymbol{k}},\omega )$ 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 ₀ in the simulations points along $\hat{z}$, 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 (the cadence at which the data are saved for analysis) was at least 2 times larger than the frequency of the fastest wave, and the total time of acquisition was longer than the period of the slowest wave in the system. It is worth mentioning that spatiotemporal spectra have been used before in numerical simulations and experiments on rotating turbulence (Clark di Leoni et al. 2014), stratified turbulence (Clark di Leoni & Mininni 2015), quantum turbulence (Clark di Leoni et al. 2015b), and incompressible MHD turbulence simulations (Meyrand et al. 2015, 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, 2010). Quantifying the relative presence of waves and/or nonlinear structures is the main outcome expected from the spatiotemporal analysis. In particular, we use the spatiotemporal spectra to search for the presence of waves and quantify the energy found near their respective dispersion relations in linear and nonlinear MHD turbulence.

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 nonrelativistic equation of motion:

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{r}}}{{dt}}={\boldsymbol{v}},\quad \displaystyle \frac{d{\boldsymbol{v}}}{{dt}}=\alpha \left({\boldsymbol{E}}+{\boldsymbol{v}}\times {\boldsymbol{B}}\right),\end{eqnarray} \tag{ 15 }$$

where the electric field E is obtained from Ohm's law, and its dimensionless (using a characteristic electric field E₀ = v₀ B₀/c) expression is

$$\begin{eqnarray}&&{\boldsymbol{E}}=\displaystyle \frac{{\boldsymbol{J}}}{{R}_{m}}-{\boldsymbol{u}}\times {\boldsymbol{B}}.\end{eqnarray} \tag{ 16 }$$

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; Balzarini et al. 2022; Pugliese & Dmitruk 2022). As such, density fluctuations will not directly affect particle motion through Equation (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₀) for particles with velocity v₀ 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₀, we can define an average gyroperiod τ_p = 2π/ω_g = 2π/α B₀. In the case of protons, we can relate α to the MHD field parameter through

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{L}_{0}}{{d}_{i}},\end{eqnarray} \tag{ 17 }$$

where ${d}_{i}={m}_{p}c/\sqrt{4\pi {\rho }_{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).

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 ${{ \mathcal 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 statistical tool we use is the standard deviation of the volumes ${\sigma }_{{ \mathcal V }}$, which should coincide with σ_RPP for a uniform distribution and increases as the preferential concentration increases.

In the aforementioned work it was found that particles accumulated in regions where ∇_⊥ · u _⊥ = ∂_x u_x + ∂_y u_y < 0. Through Equation (16) this also implies (∇ × E )_z < 0 and thus a clockwise rotating electric field. Therefore, in these regions the electric force and particle gyration are aligned and there is positive energization after a whole gyroperiod, known as betatron acceleration. This, in combination with the trapping effect of ∇_⊥ · u _⊥ < 0, makes this mechanism very effective at energizing particles.

The nonlinear simulations are performed by numerically solving Equations (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³ = 512³ 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. Here, ${k}_{d}={({\epsilon }_{d}/{\nu }^{3})}^{1/4}$ is the Kolmogorov dissipation wavenumber (which defines the dissipation scale l_d = 2π/k_d ), with ν = μ/ρ₀ the kinetic viscosity and _d the energy dissipation rate.

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 scale and timescale are defined individually for each simulation. For the characteristic length scale L₀, we use the energy-containing scale L₀ = 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}_{\mathrm{rms}}/\sqrt{4\pi {\rho }_{0}}$. For the timescales, we alternatively use the large eddy turnover time t₀ = L₀/v₀ and the particle gyroperiod τ_p = 2π/α B₀, 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₀/b_rms) of 1:9, similar to values reported in the solar wind (Hadid et al. 2017; Andrés et al. 2022).

The linear simulations are performed in the same way as the nonlinear ones but we cancel all the nonlinear terms in the MHD equations, as shown in Equations (9)–(11). In the absence of nonlinear terms, there is no energy cascade from the injection scale L₀ to smaller scales. Therefore, modes with wavenumber k outside the forced shell would quickly vanish in the presence of the dissipation terms in Equations (10) and (11). To counter this, we evolve the fields in Equations (9)–(11) without forcing and without dissipation (i.e., R_e , R_m → ∞), thus ensuring a constant energy spectrum. 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 Figure 2. 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 ${\psi }_{{\boldsymbol{k}}}\mapsto {e}^{i{\phi }_{{\boldsymbol{k}}}}{\psi }_{{\boldsymbol{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). This is because phase synchronization is a consequence of nonlinear evolution and plays an essential role in the generation of coherent structures (Kuramoto 1984). In Figure 1 are shown perpendicular slices of J_z for the initial conditions of both the L and LR runs. The L run displays current sheets as coherent structures, which are destroyed for the LR run after phase randomization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 rescaling all the Fourier coefficients, thus preserving any power-law behavior.

Initially, particles were uniformly distributed in the box with a Gaussian velocity distribution function, with an rms value of $\langle {v}_{i}^{2}{\rangle }^{1/2}\approx 1.2{v}_{0}$. As shown in Table 1, this value lies between the Alfvén velocity v₀ 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₀/l_d for all simulations. The numerical integration of Equation (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 .

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 perpendicular and parallel components to the magnetic guide field in Figures 3(a) and (b), respectively. As is known for protons interacting with MHD fields, energization is much higher in the perpendicular direction (represented here by the x component) than in the parallel direction (represented here by the z direction). For both directions, it is clear that the linear case L has significantly smaller energization than the full nonlinear simulation NL1. However, for the phase-randomizing (LR) run, the energization increases greatly and surpasses that of the 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, giving values similar to the linear L run. Furthermore, in the perpendicular case shown in Figure 3(a) all the linear runs seem to display a diffusive or subdiffusive behavior (in momentum space) $\langle {v}_{x}^{2}\rangle \sim {t}^{\alpha }$ 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 that the interaction is localized in time.

Zoom In Zoom Out Reset image size

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 Equations (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 Figure 4(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 Figure 4(b) we show the autocorrelation length ℓ_c as a function of time, where we note that in the NL1 run ℓ_c ≈ 18l_d holds 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 Section 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 Section 2.5. To further check this claim, in Figure 5 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 ${\sigma }_{{ \mathcal V }}$ could mean that this mechanism still survives in some small-scale regions, but given the diffusive energization of Figure 3 we could consider this effect secondary at best.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We turn now to the comparison of the linear energization runs amongst themselves. Having discarded interaction with structures 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 probability density function (PDF) of particle velocities at the end of each simulation. The clear geometrical difference of each PDF shown in Figure 6 suggests that each mechanism is fundamentally different, as the data are normalized to the unit variance. This normalization removes the blatant anisotropy in energization (as shown in Figure 3) and allows us to concentrate on statistical properties. 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 a constant 1-norm (i.e., ∣v_x ∣ + ∣v_z ∣ = const.), suggesting that simultaneously high perpendicular and parallel energization is unlikely, as expected from this structural energization (high-v_z particles are difficult to trap). For the L run, the shapes are clearly circular, corresponding to a constant 2-norm (${v}_{x}^{2}+{v}_{z}^{2}\,=\,$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 a constant ∞-norm ($\max \{{v}_{x},{v}_{y}\}\,=\,$const.).

Zoom In Zoom Out Reset image size

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 for LR and L in Figure 6, 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 the two 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 spatiotemporal spectra were used (e.g., Clark di Leoni et al. 2015a; Andrés et al. 2017). Therefore, we resolved the waves in time and space by choosing a very high cadence at which to store the magnetic field, in particular we used dt = 2.5 × 10⁻⁴ as the temporal sampling rate. It is worth mentioning that we assumed that the energy concentrated around the linear dispersion relation can be explained by linear and weak turbulence theories (Chandran 2005, 2008), while any spread away from the dispersion relation is a sign of strong turbulence that requires fully nonlinear theories to be understood. Figure 7 shows the spatiotemporal spectrum of the parallel magnetic field E_zz (k_x = 0, k_y , k_z = 0, ω) for the NL1 and L runs (the other linear runs have similar spectra to the L run), i.e., we separate the linear case from the nonlinear one. The dispersion relations for the fast and slow magnetosonic waves given by Equation (13) are shown as green solid and green dashed–dotted lines, respectively. Also, we include the particle gyrofrequency ω_g , which lies in the wavenumber range 40 < ∣ k ∣ < 80, and the characteristic sweeping frequency given by ${\omega }_{\mathrm{sw}}\sim {U}_{\mathrm{rms}}\sqrt{{k}_{\perp }^{2}+{k}_{\parallel }^{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 or nonpropagating modes) for all wavenumbers. We also noted a small amount 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.

Zoom In Zoom Out Reset image size

Furthermore, the absence of sweeping in the linear case is also useful for understanding the low energization in the L run. As shown in Figure 4, 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.

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 Figure 3), but not quantitatively. The same occurs for the correlation length and volume deviation ${\sigma }_{{ \mathcal V }}$, as summarized in Table 2 (energization and ${\sigma }_{{ \mathcal V }}$ are calculated at the end of the simulation). We see that both energization and correlation length decrease along with τ_c , while ${\sigma }_{{ \mathcal 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) prevents structures from being stable, for lack of either size or intensity. Second, we could conclude that particles tend to accumulate similarly in all cases, but their ability to exploit energization mechanism diminishes for lower τ_c .

Table 2. Particle-related Quantities in Nonlinear Simulations

Run	τc /τp	ℓc /ld	$\langle {\rm{\Delta }}{v}_{x}^{2}\rangle /{v}_{0}^{2}$	${\sigma }_{{ \mathcal V }}/{\sigma }_{{\rm{RPP}}}$	$\langle { \mathcal S }\rangle /{\tau }_{p}$
NL1	1.146 × 101	18.0	24.3	1.87	2.23
NL2	1.146 × 100	15.2	20.3	1.98	2.13
NL3	2.865 × 10−1	14.4	15.6	1.70	1.85

Note. Quantities shown are forcing correlation time, field correlation length, mean value of the perpendicular energization, volume deviation, and clustered time of the particles for nonlinear simulations. While correlation length, energization, and clustered time decrease with τ_c , volume deviation remains mostly constant.

Download table as:  ASCIITypeset image

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 ${ \mathcal V }$ is lower than some threshold value ${{ \mathcal V }}_{{\rm{th}}}$, which is obtained by comparing the PDF of the cell volume with that of an RPP (see references in Section 2.5 for more details). As seen in Figure 5, clustering takes some time to settle, and as such the label "clustered" may not mean much initially. Toward the end of the simulation, "clustered" is almost equivalent to "trapped inside a structure." Then, we can determine when and where particles are accumulating and therefore calculate the continuous time they spent clustered. We define a streak as an interval of (discrete) time for which the particle is clustered and its duration as ${ \mathcal S }$. In particular, we calculated the mean streak time $\langle { \mathcal S }\rangle$ 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 two 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 unclustered particles) in low- τ_c simulations, suggesting clusters become more feeble and unstable with shorter correlation time. For smaller 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 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 Figure 8. We can see that after ∼100τ_p (long enough for clustering to settle), energization becomes exponential in time.

Zoom In Zoom Out Reset image size

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. Using Equation (15), we could derive an evolution equation for the perpendicular energy ε_⊥,

$$\begin{eqnarray}&&\displaystyle \frac{d{\varepsilon }_{\perp }}{{dt}}=\alpha \left[{{\boldsymbol{E}}}_{\perp }\cdot {{\boldsymbol{v}}}_{\perp }-\left({{\boldsymbol{v}}}_{\parallel }\times {{\boldsymbol{v}}}_{\perp }\right)\cdot {{\boldsymbol{b}}}_{\perp }\right]\equiv {{ \mathcal P }}_{\perp }-{{ \mathcal P }}_{\times },\end{eqnarray} \tag{ 18 }$$

with α defined in Equation (17). The first term ${{ \mathcal P }}_{\perp }=\alpha {\boldsymbol{E}}\cdot {{\boldsymbol{v}}}_{\perp }$ corresponds to net energization while the second relates to exchange with parallel energy, such as pitch-angle scattering. Averaging over one gyroperiod,

$$\begin{eqnarray}\begin{array}{rcl}\langle {{ \mathcal P }}_{\perp }{\rangle }_{{\tau }_{p}} & = & \displaystyle \frac{1}{{\tau }_{p}}{\int }_{0}^{{\tau }_{p}}\alpha {{\boldsymbol{E}}}_{\perp }\cdot {{\boldsymbol{v}}}_{\perp }dt\approx \displaystyle \frac{\alpha }{{\tau }_{p}}{\oint }_{{ \mathcal C }}{{\boldsymbol{E}}}_{\perp }\cdot {dl}\\ & = & -\displaystyle \frac{\alpha }{{\tau }_{p}}{\iint }_{\left.{ \mathcal S }({ \mathcal C }\right)}{\boldsymbol{\nabla }}{E}_{\perp }\cdot {dS},\end{array}\end{eqnarray} \tag{ 19 }$$

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 Equation (16) as E _⊥ ≈ – u _⊥ × B ₀ and therefore the parallel component of its curl is ${({\boldsymbol{\nabla }}{E}_{\perp })}_{\parallel }\approx {B}_{0}{{\boldsymbol{\nabla }}}_{\perp }\cdot {{\boldsymbol{u}}}_{\perp }$. Expanding the integral at the lowest order in the gyroradius R_g , we achieve

$$\begin{eqnarray}&&\langle {{ \mathcal P }}_{\perp }{\rangle }_{{\tau }_{p}}\approx -\displaystyle \frac{{\alpha }^{2}{B}_{0}^{2}}{2\pi }{\iint }_{\left.{ \mathcal S }({ \mathcal C }\right)}{{\boldsymbol{\nabla }}}_{\perp }\cdot {{\boldsymbol{u}}}_{\perp }dS\approx -\displaystyle \frac{{\alpha }^{2}{B}_{0}^{2}}{2\pi }\pi {R}_{g}^{2}{{\boldsymbol{\nabla }}}_{\perp }\cdot {{\boldsymbol{u}}}_{\perp }.\end{eqnarray} \tag{ 20 }$$

Using R_g = ∣ v _⊥∣/α B₀ and recalling the definition of ε_⊥, we can substitute in Equation (18) to obtain

$$\begin{eqnarray}&&\displaystyle \frac{d{\varepsilon }_{\perp }}{{dt}}\approx \lambda {\varepsilon }_{\perp },\quad \lambda =-{{\boldsymbol{\nabla }}}_{\perp }\cdot {{\boldsymbol{u}}}_{\perp },\end{eqnarray} \tag{ 21 }$$

which for approximately constant ∇_⊥ · u _⊥ < 0 predicts an exponential increase for ε_⊥. As pointed out in Section 2.5, this is produced by the alignment of electric field and particle velocity and could be understood as betatron-like acceleration (Swann 1933) following the plasma flow, which trapped particles naturally do.

The slopes in Figure 8 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 $\langle { \mathcal S }\rangle$. Particles can leave structures in basically three 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₀ 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}_{\perp }^{2}\sim {(\alpha {B}_{0}{{\ell }}_{c})}^{2}\sim 1600{v}_{0}^{2}$, which is achieved by practically no particle (less than 1 particle in 10⁴) 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.

For the study of the relevance of waves in the system for different τ_c , we quantitatively analyze the spatiotemporal magnetic spectra. Figure 9 shows the spatiotemporal 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 show the dispersion relation for the Alfvén, fast, and slow magnetosonic waves as blue dashed, green dashed–dotted, and orange solid lines, respectively. We also added the particle gyrofrequency and the sweeping frequency as red dashed and blue solid lines, respectively. We observed that the energy is mainly located around the slow branch and, to a 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 toward the sweeping frequency. This result is indeed compatible with the behavior of ℓ_c /l_d , as the sweeping energy is mostly related to the nonlinear structures (see Table 2).

Zoom In Zoom Out Reset image size

In order to analyze the magnetosonic modes, we studied the spatiotemporal spectrum of a parallel magnetic field fluctuation E_zz (k_x = 0, k_y , k_z = 15). Figure 10 shows the same trend that we observed in Figure 9. In fact, the energy around the fast branch decreases noticeably as we increase τ_c . It is worth mentioning that the slow branch is completely immersed in the sweeping region, therefore we cannot come to a conclusion about the increase (or decrease) in the magnetic energy around the slow mode. However, in the case of NL1, the energy is more concentrated around ω = 0 than for smaller values of the correlation time, and the sweeping frequency correctly delimits this energy.

Zoom In Zoom Out Reset image size

For a more precise study, we used an integration method to quantify the amount of energy near the different wave modes in each spatiotemporal 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

$$\begin{eqnarray}&&F({k}_{z})=\displaystyle \frac{\left.{E}_{{xx}}({k}_{x}=0,{k}_{y}=0,{k}_{z},\omega ={\omega }_{{\rm{A}},{\rm{F}},{\rm{S}}}\right)}{\left.{{\rm{\Sigma }}}_{j}{E}_{{xx}}({k}_{x}=0,{k}_{y}=0,{k}_{z},{\omega }_{j}\right)},\end{eqnarray} \tag{ 22 }$$

with ω_A,F,S the frequencies that satisfy a certain dispersion relation (and the E_xx (k_z ) component is an illustrative example only). Figure 11 shows the energy around (a) Alfvén and (b) slow magnetosonic waves for runs NL1, NL2, and NL3. In particular, we observe 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 in the rest of the simulations as a result of the energy moving toward the sweeping frequency. In the case of Alfvén modes, most of the energy is concentrated around low wavenumbers and it 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 12 shows (a) the energy near fast magnetosonic waves for runs NL1, NL2, and NL3 and (b) the difference between NL3 and the first two nonlinear 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%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size