© The Authors 2024

1. Introduction

Plasma turbulence in the presence of a sufficiently strong background magnetic field B₀ exhibits a spectral anisotropy: fluctuating energy is concentrated in modes with strongly oblique wave-vectors with respect to B₀ (Montgomery & Turner 1981; Shebalin et al. 1983; Ghosh & Goldstein 1997). This anisotropy is not well understood; we have phenomenological approaches such as those based on geometric decompositions, for example two-dimensional (2D) plus slab, or conjectures like critical balance (Matthaeus et al. 1990; Goldreich & Sridhar 1995; Oughton et al. 2015; Oughton & Matthaeus 2020). To avoid the need of such ad hoc assumptions we can characterise turbulence and its anisotropy using the Kármán–Howarth–Monin (KHM) equation starting from first principles. It was derived in the incompressible, constant density approximation, first for hydrodynamic turbulence (de Kármán & Howarth 1938; Monin & Yaglom 1975), and then extended to magneto-hydrodynamic (MHD) turbulence (Politano & Pouquet 1998a,b) and Hall MHD (Galtier 2008; Hellinger et al. 2018; Ferrand et al. 2019). This equation quantifies different processes. In a decaying (non-driven) system, this equation gives the rate of losses (decay) of energy at large scales (in the energy-containing range), the cross-scale energy transfer (cascade), as well as the dissipation. This equation is formulated in the separation space of increments, which gives a natural possibility for investigation of turbulence anisotropy in the separation space looking at increments parallel and perpendicular to the background magnetic field B₀.

Results of the KHM analysis of direct simulation of MHD turbulence (Verdini et al. 2015) show a clear anisotropy of structure functions in the separation space (l_⊥, l_∥). This work also highlights the somewhat paradoxical properties of turbulence in magnetised plasmas. The inertial range can be formally defined as the region where the cascade rate K is equal to the dissipation rate Q. The simulation results of Verdini et al. (2015) clearly exhibit a tendency to have K(l_⊥, l_∥) = Q within a range of scales that is anisotropically distributed in the (l_⊥, l_∥) space. In the isotropic case K just expresses the cascade rate, while in the anisotropic case this quantity characterises the cascade rate but does not have any information about the direction of energy cascade in the (l_⊥, l_∥) space.

To describe the cross-scale energy transfer in turbulence spectral or spatial filtering (Eyink & Aluie 2009; Aluie 2011; Grete et al. 2017) techniques are also used. In this case, however, filters are usually assumed to be isotropic; when anisotropic spatial filtering (coarse graining) is used the cascade rate exhibits an anisotropy (Manzini et al. 2022) that is similar to the KHM results.

For weakly collisional plasmas (such as those in the solar wind), it is necessary to include the effects of plasma compressibility (Hellinger et al. 2021a) and treat the pressure as a tensor (Hellinger et al. 2022). The pressure–strain interaction, which couples the magnetic+kinetic energy to the internal energy (Del Sarto et al. 2016), may work as an effective dissipation channel, as seen in direct kinetic numerical simulations (Yang et al. 2017; Matthaeus et al. 2020; Hellinger et al. 2022; Arró et al. 2022).

Numerical simulations and in situ observations suggest that the pressure–strain channel is likely responsible for correlations between particle velocity-field gradients and temperatures (Franci et al. 2016; Parashar & Matthaeus 2016; Pezzi et al. 2021; Yordanova et al. 2021), but its properties are not well understood. Many different physical processes contribute to this interaction (Yang et al. 2017; Cassak & Barbhuiya 2022), including those connected with plasma compression and expansion. It is not clear at what spatial scales the pressure–strain effect acts. Electron and ion energisations are expected to work at very different scales (Franci et al. 2022), but comparable scales are observed in some full-particle simulations (Yang et al. 2022). It is also unclear what the role of the turbulence anisotropy is (most of the present simulations are two-dimensional) and how it compares to other processes such as Hall coupling (Papini et al. 2019, 2021). In this paper we extended the two-dimensional work of Hellinger et al. (2022), we analysed a three-dimensional hybrid simulation of decaying plasma turbulence using the KHM equation, and we studied the anisotropy of different turbulent processes.

2. The Kármán–Howarth–Monin equation

2.1. Extended Hall magnetohydrodynamics

We investigated a system governed by the following extended compressible inviscid Hall MHD equations for the plasma density ρ, the plasma mean velocity u, and the magnetic field B:

$$\begin{array}{cc}& {\partial }_t\rho +(\mathit{u}·\mathbf{\nabla })\rho =-\rho \mathbf{\nabla }·\mathit{u},\hfill \\ \hfill & \rho {\partial }_t\mathit{u}+\rho (\mathit{u}·\mathbf{\nabla })\mathit{u}=\mathit{J}\times \mathit{B}-\mathbf{\nabla }·\mathbf{P},\hfill \\ \hfill & {\partial }_t\mathit{B}=\mathbf{\nabla }\times [(\mathit{u}-\mathit{j})\times \mathit{B}]+\eta {\mathrm{\nabla }}^2\mathit{B}.\hfill \end{array}$$(1)

Here P is the plasma pressure tensor, η is the electric resistivity, J is the electric current density, j is the electric current density in velocity units, and j = J/ρ_c = u − u_e (ρ_c and u_e being the ion charge density and the electron velocity, respectively). We assume SI units except for the magnetic permeability μ₀, which is set to one (SI results can be obtained by the rescaling $\mathit{B}\to \mathit{B}{\mu }_0^{-1/2}$). Equation (1) (except the resistive term) can be derived taking moments of the Vlasov equation for protons and electrons, and assuming massless electrons. We added a resistive dissipation, as used in the hybrid code to avoid the accumulation of energy on small scales; this process substitutes the electron dissipation. The plasma pressure tensor contains contributions from protons and electrons but for electrons we further adopted a scalar pressure (an approximation used in the hybrid code).

For the sum of kinetic (E_kin = ⟨ρ|u|²/2⟩) and magnetic (E_mag = ⟨|B|²/2⟩) energies (averaged over a closed volume) we obtained this budget equation,

$$\begin{array}{c}\hfill {\partial }_t(E_{\mathrm{kin}}+E_{\mathrm{mag}})=⟨\mathbf{P}:\mathbf{\nabla }\mathit{u}⟩-{\eta ⟨|\mathit{J}|}^2⟩,\end{array}$$(2)

where the colon (:) operator denotes a double contraction of two tensors (i.e., for two tensors T and V, T : V = ∑_ijT_ijV_ij, T_ij and V_ij being their components) and the angle brackets ⟨ • ⟩ denote the averaging. The kinetic and magnetic energy is transferred to the internal energy by the resistive dissipation rate

$$\begin{array}{c}\hfill Q_{\eta }={\eta ⟨|\mathit{J}|}^2⟩.\end{array}$$(3)

The exchanges between the kinetic+magnetic energy and the internal energy also proceed through the pressure–strain coupling. We defined the pressure–strain rate as

$$\begin{array}{c}\hfill \psi =-⟨\mathbf{P}:\mathbf{\nabla }\mathit{u}⟩.\end{array}$$(4)

Since the pressure–strain may behave as an effective dissipation rate, we defined the total (effective) dissipation rate as a sum of the two:

$$\begin{array}{c}\hfill Q=Q_{\eta }+\psi .\end{array}$$(5)

We did not consider an explicit viscosity since it is not present in the hybrid approximation.

2.2. The KHM equation

For the formulation of the KHM equation in terms of structure functions for compressible Hall MHD, we used the density-weighted velocity field w = ρ^1/2u (Kida & Orszag 1990) to take into account a variable density. We characterised the spatial scale decomposition of the kinetic and magnetic energies (and their sum) using the second-order structure functions

$$\begin{array}{c}\hfill S_w=\frac{1}{4}\langle {|\delta \mathit{w}|}^2\rangle ,S_B=\frac{1}{4}\langle {|\delta \mathit{B}|}^2\rangle ,\mathrm{and}S=S_w+S_B,\end{array}$$

where the deltas denote increments of the corresponding quantities (e.g., δw = w(x + l)−w(x), x being the position and l the separation vector).

For the second-order structure function S representing the sum of the kinetic and magnetic energies, we derived the KHM equation in the form

$$\begin{array}{c}\hfill {\partial }_{t}S=K_{\mathrm{MHD}}+K_{\mathrm{Hall}}-\mathrm{\Psi }-D,\end{array}$$(6)

where K_MHD and K_Hall are the MHD and Hall cascade rates, respectively; Ψ represents the pressure–strain effect; and D accounts for the effects of dissipation and heating. The cascade rates K_MHD and K_Hall (as well as Ψ and D) can be expressed in a few different equivalent forms. Here we used a form that uses increments of second-order quantities (equivalent to the usual form involving divergence of a third-order structure function vector; cf. Hellinger et al. 2021a):

$$\begin{array}{cc}\hfill K_{\mathrm{MHD}}& =-\frac{1}{2}\langle \delta \mathit{w}·\delta [(\mathit{u}·\mathbf{\nabla })\mathit{w}+\frac{1}{2}\mathit{w}(\mathbf{\nabla }·\mathit{u})]\rangle \hfill \\ \hfill & -\frac{1}{2}\langle \delta \mathit{w}·\delta \left(\frac{\mathit{J}\times \mathit{B}}{\sqrt{\rho }}\right)\rangle +\frac{1}{2}\langle \delta \mathit{J}·\delta (\mathit{u}\times \mathit{B})\rangle ,\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill K_{\mathrm{Hall}}& =-\frac{1}{2}\langle \delta \mathit{J}·\delta (\mathit{j}\times \mathit{B})\rangle ,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \mathrm{\Psi }& =-\frac{1}{2}\langle \delta \mathit{w}·\delta \left(\frac{\mathbf{\nabla }·\mathbf{P}}{\sqrt{\rho }}\right)\rangle ,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill D& =\frac{1}{2}\eta \langle {|\delta \mathit{J}|}^2\rangle =Q_{\eta }-\frac{1}{2}\eta {\mathrm{\nabla }}^2{S}_B.\hfill \end{array}$$(10)

3. Simulation results

Here we analysed results of three-dimensional (3D) hybrid simulations of decaying plasma turbulence. In the hybrid approximation, ions are described by a particle-in-cell model, whereas electrons are a massless charge neutralising fluid (Matthews 1994; Franci et al. 2018a,b)¹ The system is numerically integrated in a 3D domain (x, y, z) with 512³ grid points and the spatial resolution Δx = Δy = Δz = d_i/8. Here d_i denotes the ion inertial length. Protons are initially isotropic with β_i = 0.5, β_i being the ion beta (i.e., the ratio of the ion to magnetic pressures). In order to reduce the noise, a Gaussian smoothing on 3 × 3 points is used on the proton density and velocity in the code. A uniform ambient magnetic field B₀, directed along z is present, whereas neutralising electrons are assumed to be isothermal with β_e = β_i; here β_e denotes the electron beta, the ratio of the electron and magnetic pressures. The system was initialised with an isotropic 3D spectrum of modes with random phases, linear Alfvén polarisation (δB ⊥ B₀), and vanishing correlation between magnetic field and velocity fluctuations. These modes were in the range $k\le 0.22d_{\mathrm{i}}^{-1}$ with rms fluctuations δB = 0.25. The time step is $\mathrm{\Delta }t=0.005{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$ for particle integration (the magnetic field is advanced with a smaller time step Δt_B = Δt/20), the number of particles per cell N_ppc = 4096, and a small resistivity $\eta ={10}^{-3}{\mu }_0{v}_{\text{A}}^2/{\mathrm{\Omega }}_{\text{i}}$ is used to avoid energy accumulation at the smallest scales; we note that no explicit viscosity is present in the hybrid model. Here Ω_i denotes the ion cyclotron frequency, μ₀ is the magnetic permeability of the vacuum, and v_A stands for the Alfvén velocity. We let the system evolve until turbulence becomes quasi-stationary.

In the simulation, the proton kinetic energy E_kin and the magnetic energy E_mag oscillate with opposite phases, thus suggesting energy exchanges. However, overall these two quantities decrease. On the other hand, the proton internal energy, E_int = 3⟨p⟩/2 increases (here p is the scalar pressure, p = tr(P), tr being the trace). This energisation works throughout the simulation even before turbulence is well developed. The total energy E_tot = E_kin + E_int + E_mag slowly decreases owing to the resistive dissipation since electrons are assumed to be massless and isothermal. This behaviour is shown on Fig. 1a, which displays the evolution of relative changes of the kinetic, magnetic, and internal energies ΔE_{kin, mag, int, tot} given by ΔE(t) = [E(t)−E(0)]/E_tot(0). The proton energisation is driven by the pressure–strain coupling that plays the role of an effective dissipation channel. Figure 1b shows the resistive dissipation rate Q_η as a function of time as well as the two contributions to the pressure–strain rate ψ, the compressive contribution (Yang et al. 2017)

$$\begin{array}{c}\hfill {\psi }_{\mathrm{c}}=-⟨p\mathbf{\nabla }·\mathit{u}⟩,\end{array}$$(11)

Fig. 1. Evolution of different quantities as a function of time: (a) Relative changes in the kinetic energy ΔEkin (dashed line), magnetic energy ΔEmag (solid line), internal energy ΔEint (dash-dotted line), and total energy ΔEtot (dotted line); (b) resistive dissipation rate Qη (dashed line), incompressive pressure–strain effective dissipation rate ψi (solid line), and its compressive counterpart ψc (dotted line). These rates are normalised to the total effective dissipation rate Q = Qη + ψ (averaged over 200 ≤ tΩi ≤ 300).

and the incompressive contribution ψ_i = ψ − ψ_c. These terms are normalised to the total effective dissipation rate $Q=7.15\times {10}^{-5}{m}_{\text{p}}{v}_{\text{A}}^2{\mathrm{\Omega }}_{\text{i}}$ obtained by averaging ψ + Q_η over the time interval 200 ≤ tΩ_i ≤ 300. The resistive dissipation rate Q_η starts at zero, slowly increases, and becomes quasi-stationary for tΩ_i ≳ 200. The effective dissipation rate ψ_i is initially dominant, and strongly and rapidly oscillates in time for tΩ_i ≲ 50; later on its oscillations become weaker and slower. At later times Q_η and ψ_i become comparable to ψ; during the interval 200 ≤ tΩ_i ≤ 300 the averaged ψ_i ≃ 0.47Q, whereas Q_η ≃ 0.44Q. The compressive component of the pressure–strain interaction ψ_c starts at zero (owing to the incompressive initialisation) and reaches values comparable to Q. At later times it oscillates around zero, so that its averaged value during the time interval 200 ≤ tΩ_i ≤ 300 is ψ_c ≃ 0.08Q.

The non-linear coupling between the magnetic and plasma velocity fields leads to a turbulent cascade. The value of Q_η becomes quasi-stationary for t ≳ 170Ω_i and we assume that the system reaches well-developed turbulence (Mininni & Pouquet 2009; Servidio et al. 2015).

4. Isotropic view

The turbulent cascade leads to a spread of fluctuating magnetic and kinetic energies over a wide range of scales. During the last quasi-stationary phase their spectral properties are only weakly varying. An example is given in Fig. 2 that shows omnidirectional power spectral densities of the magnetic field B, P_B (in red), and the proton velocity field w = ρ^1/2u, P_w (in blue), and their sum P (in black) as a function of k. The magnetic spectrum is roughly Kolmogorov-like (∝k^−5/3) at large scales with a break at around kd_i ∼ 3 where it steepens. The proton kinetic spectrum has a somewhat steeper slope at large scales compared to P_B with a break at around kd_i ∼ 1 where it strongly steepens. For kd_i ≳ 5 the kinetic spectrum is affected by the numerical noise connected with the finite number of particles per cell. The magnetic spectrum is also influenced by the noise and this effect is important for kd_i ≳ 10.

Fig. 2. Omnidirectional power spectral densities of the magnetic field B, PB (red), and the compensated proton velocity field w, Pw (blue), and their sum P (black) as a function of k normalised to di at $t=250{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$. The dotted line shows a spectrum ∝k−5/3 for comparison.

We can test turbulence properties in more detail using the KHM equation. While the background magnetic field leads to an important spectral anisotropy (as we discuss later), for the assessment of the overall turbulence behaviour we start with isotropised quantities. Before we present the analysis of the simulation results via the KHM equation, we start with a presentation of the theoretical expectations. We assume a simplified model described with

$$\begin{array}{c}\hfill {\partial }_{t}S=K-D,\end{array}$$(12)

where ∂_tS represents the decay of the energy, K is the cascade rate, and D stands for the dissipation. The KHM equation has cumulative properties (Hellinger et al. 2021a,b, see also Appendix A), the two terms, ∂_tS(l) and D(l), behave as a spatial low-pass filter and correspond to the decay and dissipation rates on scales smaller than l. In contrast, the term K describes the cross-scale energy transfer (cascade) rate. Figure 3 shows the expected behaviour of the KHM terms as a function of l. The value of D starts at zero and increases towards Q. Above the dissipation length l_diss, where dissipation is not active, D is then constant: D = Q. In the ideal decaying case the system is stationary in the dissipation range and in the inertial range, and ∂_tS is zero there. Only on large scales l > l_outer does this term become negative, ∂_tS → −Q (this is the energy containing range). Finally, between the dissipation length l_diss and the outer scale l_outer the cascade rate is constant, equal to the dissipation rate K = Q (this relation is the formal definition of the inertial range).

Fig. 3. Schematic view of the theoretical expectations of a simplified system 0 = −∂tS + K − D. The decay term −∂tS (blue), the cascade rate K (green), and the dissipation rate −D (red) are given as a function of l and are normalised to the dissipation rate Q.

The KHM equation constitutes a cross-scale energy conservation equation, and its application represents a strong validity check: one cannot expect that the KHM equation holds exactly in a numerical code. We defined the validity test O,

$$\begin{array}{c}\hfill O=-{\partial }_{t}S+K_{\mathrm{MHD}}+K_{\mathrm{Hall}}-\mathrm{\Psi }-D,\end{array}$$(13)

which measures the error of the code owing to numerical and other issues. Figure 4 displays the validity test and the different contributing terms (−∂_tS, K_MHD, K_Hall, −Ψ, −D) normalised to the total effective dissipation rate Q as a function of the scale separation l = |l| (averaged over the quasi-stationary time period 200 ≤ tΩ_i ≤ 300). Figure 4a shows that the KHM Eq. (6) is relatively well satisfied. The simulated system exhibits the expected behaviour, but is far from the ideal schema of Fig. 3. In this case there is no clear separation between the dissipation and the energy containing ranges. The energy decay ∂_tS is negative and important over a wide range of scales; it becomes negligible somewhat below l = d_i. At the same time, the dissipation processes are active to larger scales. Consequently, the total cascade K_MHD + K_Hall rate only reaches a maximum value of about one-half the dissipation rate Q. Figure 4b shows that the resistive dissipation is active on small scales l ≲ 2d_i, but that the pressure–strain interaction is active at relatively large scales, l ≲ 20d_i, largely extending the dissipation range. On intermediate scales the dominant process is the MHD cascade; K_MHD is negative at large scales, which indicates an inverse cascade. On small scales part of the cascade continues via the Hall term, K_Hall attains a maximum value of about 15% of Q.

Fig. 4. Isotropised KHM analysis averaged results over a period of pressure–strain oscillations (200 ≤ tΩi ≤ 300): (a) Validity test O (black line) as a function of l along with combined dissipative terms −(D + Ψ) (red), combined cascade rates KMHD + KHall (green), and the decay rate −∂tS (blue) in a format similar to Fig. 3. (b) Different separate contributing terms to the KHM equation: the decay rate −∂tS (blue), the MHD cascade rate KMHD (green), the Hall cascade rate KHall (orange), the resistive dissipation rate −D (red), and the pressure–strain (effective dissipation) rate −Ψ (magenta). The dashed magenta curve displays the compressive component of −Ψ. All the quantities are normalised to the time-averaged effective total dissipation rate Q.

Compressive effects contribute to the pressure–strain coupling and can be described as

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{c}}=-\frac{1}{2}\langle \delta \mathit{w}·\delta \left(\frac{\mathbf{\nabla }p}{\sqrt{\rho }}\right)\rangle \end{array}$$(14)

(see Eq. (9)). The magenta dashed line in Fig. 4 shows the Ψ_c contribution to Ψ as in the case of the pressure-dilatation term ψ_c, Ψ_c is weak when averaged over 200 ≤ tΩ_i ≤ 300.

The KHM equation is valid during the whole simulation (it is relatively well satisfied |O|/Q ≲ 10%), so we can investigate how turbulence begins, evolves, and varies. Figure 5 shows the evolution of the isotropised KHM results, the different contributing terms as a function of time, and the separation scale l = |l| normalised to the total effective dissipation rate Q (averaged over the last third of the simulation 200 ≤ tΩ_i ≤ 300).

Fig. 5. Evolution of isotropised KHM results. Shown are the different KHM terms as a function of time t and l: (a) decay rate ∂tS, (b) MHD cascade rate KMHD, (c) Hall cascade rate KHall, (d) resistive dissipation rate D, (e) compressive pressure–strain rate Ψc, and (f) incompressive pressure–strain rate Ψi. All the quantities are normalised to the effective total dissipation rate Q.

Figure 5 demonstrates that the features of well-developed turbulence seen in Fig. 4 only appear at later times. During the turbulence onset (for about tΩ_i ≲ 200) ∂_tS is positive on intermediate scales. This means that the kinetic+magnetic energy increases in time on these scales owing to the MHD non-linear term (since ∂_tS ∼ K_MHD). At later times, ∂_tS becomes negative, monotonically decreasing with l as expected for the cumulative low-pass filter property. The range where ∂_tS is positive shifts with time to smaller scales as these are produced. As the energy transfer towards smaller scales continues, the resistive dissipation gradually appears. At later times we recover the cumulative properties of the resistive term D, increasing from 0 to about a constant value for l ≳ 2d_i. The cross-scale energy transfer (cascade) rate K_MHD remains negative on large scales (which is a possible signature of an inverse cascade), but it is positive on intermediate scales indicating an energy transfer (cascade) from large to small scales. The range of positive K_MHD slowly shifts to smaller scales over time as these are produced. For 100 ≲ tΩ_i ≲ 200 positive values of K_Hall gradually appear as the Hall physics sets in.

The evolution of the pressure–strain interaction is more complicated. As in the case of pressure–strain rate, We split Ψ into its compressive part Ψ_c and the remaining Ψ_i = Ψ − Ψ_c, incompressive part (corresponding to the Pi-D term of Yang et al. 2017). The compressive contribution Ψ_c exhibits oscillations between about −0.5Q and 0.5Q within a range of scales that extends to smaller scales as the cascade develops. Oscillations between positive and negative values indicate reversible exchanges between the magnetic+kinetic energy and the internal energy; as can be seen in Fig. 4, the average value of Ψ_c during tΩ_i ≳ 200 is much weaker than its amplitude of oscillation. In contrast, the incompressive part Ψ_i remains mostly positive; however, small negative values of both Ψ_i and Ψ_c for tΩ_i ≲ 50 suggest some weak opposite transfer on small scales. After a transient period (for tΩ_i ≳ 30), Ψ_i monotonically increases with l (as expected for the cumulative low-pass filter) similarly to D; it works as an effective dissipation mechanism. The range where Ψ_i is active (i.e., where it strongly increases with l) shifts with time from large to intermediate scales. The system is quasi-stationary for tΩ_i ≳ 200; the properties seen in Fig. 4 do not qualitatively change, but only exhibit some oscillations due to the compressive pressure–strain component Ψ_c.

We started with the spectral anisotropy of turbulent fluctuations. Figure 6 displays the (gyrotropised) power spectral densities of the magnetic energy P_B and the kinetic energy P_w as a function of k_⊥ and k_∥. Figure 6 shows that the power spectra are strongly anisotropic in the (k_⊥, k_∥) plane, indicating a preferential cascade direction that is strongly oblique with respect to the ambient magnetic field B₀. Similarly to the isotropised spectra, P_w decreases with k_⊥ faster than P_B. Their spectral anisotropies are weak on large scales (due to the initial isotropic spectrum) and increase towards small scales. The magnetic power spectrum P_B exhibits a change in behaviour for k_⊥d_i ≳ 1, a possible indication of the effect of the Hall term.

Fig. 6. Color scale plots of gyrotropised power spectral densities of (a) the magnetic energy PB, and (b) the kinetic energy Pw as a function of k⊥ and k∥ (normalised to di).

To analysed the different turbulent processes we used the full KHM equation where we gyrotropised the terms to obtain two-dimensional results as a function of l_∥ = |l ⋅ b₀| and l_⊥ = |l − (l ⋅ b₀)b₀| (where b₀ = B₀/|B₀|). Figure 7 displays the anisotropic properties of the KHM equation, iso-contours of different quantities as a function of the separation scales l_⊥ and l_∥ (averaged over the time period 200 ≤ tΩ_i ≤ 300). Figure 7a shows the second-order structure function S as a reference (along with a circle $l_{\perp }^2+l_{\Vert }^2=\text{const}.$ to emphasise the anisotropy and the effect of the log-log plot); S decreases from large to small scales and exhibits a clear anisotropy with respect to the ambient magnetic field; at large scales S is weakly anisotropic (due to the isotropic initialisation) and the anisotropy increases towards smaller scales in agreement with the spectral properties of fluctuating magnetic and kinetic energies (see Fig. 6). The different terms in the KHM equation exhibit different anisotropic properties in the (l_⊥, l_∥) plane, but their overall behaviours are similar to that of S: their anisotropy is weak on large scales and strong on small scales. Figure 7b shows the decaying term ∂_tS that is important on large scales; its main anisotropy is that |∂_tS| decreases faster with decreasing l_∥. Figure 7c displays the MHD cascade rate K_MHD; the region where K_MHD dominates shifts to larger scales and shrinks when moving from the quasi-perpendicular to quasi-parallel angles. The sign of K_MHD becomes negative at large scales, which indicates some (anisotropic) inverse cascade. Figure 7d shows the Hall cascade rate K_Hall; the Hall cascade rate exhibits a strong anisotropy, the Hall cascade (and the Hall term) becomes important at larger scales for quasi-parallel angles compared to quasi-perpendicular ones. The shrinking and shifting of the MHD cascade rate K_MHD becomes less pronounced (but does not disappear) when we consider the total cascade rate K = K_MHD + K_Hall. Figure 7e displays the pressure–strain term Ψ. It exhibits a similar anisotropy to the Hall cascade rate on small scales but becomes more isotropic on intermediate scales; the range of scales where Ψ varies is narrower at quasi-parallel scales compared to the quasi-perpendicular scales. Finally, Fig. 7f shows the resistive dissipation term; it exhibits a similar strong anisotropy to the Hall cascade rate K_Hall.

Fig. 7. Anisotropic KHM results, labelled iso-contours of different quantities as a function of the separation scales l⊥ and l∥ (averaged over the time period 200 ≤ tΩi ≤ 300): (a) Second-order structure function S, (b) Decay rate ∂tS, (c) MHD cascade rate KMHD, (d) Hall cascade rate KHall, (e) pressure–strain (effective dissipation) term Ψ, and (f) resistive dissipation term D. The dotted line in panel a displays a circle $l_{\perp }^2+l_{\Vert }^2=\text{const}.$ as a measure of anisotropy. The dotted lines in panel c denote non-positive values. The KHM quantities are normalised to the (time-averaged) effective total dissipation rate Q.

The compressive effects that lead to large oscillations of Ψ_c lead to important and anisotropic oscillations of ∂_tS, K_MHD, and Ψ mainly on large scales. However, further analysis of this phenomenon is beyond the scope of this paper.

5. Discussion

In this paper we studied decaying plasma turbulence using the 3D hybrid code. The spectral properties of the simulated system are similar to those in previous works (Markovskii & Vasquez 2011; Franci et al. 2018b). We concentrated our analysis of turbulence on the KHM equation (in the form that can be easily calculated using the fast Fourier transform Banerjee & Galtier 2017; Hellinger et al. 2021a; Montagud-Camps et al. 2022). The KHM results show how turbulence evolves from large-scale zero-cross helicity Alfvénic-like fluctuation. The non-linear coupling leads to a transfer (cascade) of energy to small scales and the resistive dissipation rate slowly increases until it becomes quasi-stationary. On the other hand, the pressure–strain rate starts to be important early in the simulation and continues to be strong (i.e., of the same order as the resistive dissipation rate) during the whole simulation and exhibits important temporal variations. These oscillations are caused by the compressive part of the pressure–strain interaction that leads to largely reversible energy exchanges between the kinetic+magnetic energy and the internal energy; analogous behaviour is observed in hydrodynamic simulations (Hellinger et al. 2021b). In contrast, the incompressive (Pi-D) part of the pressure–strain coupling works as an effective dissipation mechanism. Our results suggest that the compressibility effects tend to be negligible for a good statistical sample; in our case we took a spatio-temporal average, these effects may be important locally and may affect for instance the cascade rates estimated from observations based on limited time intervals (Marino & Sorriso-Valvo 2023).

The onset of turbulence quantified using the KHM equation is similar to the results of the 2D Hall MHD simulations (Hellinger et al. 2021a), except that the pressure–strain effective dissipation appears early in the simulation. After the system reaches a quasi-stationary state it exhibits an overall behaviour similar to that of 2D hybrid simulations: the fluctuating (kinetic+magnetic) energy decreases at large scales; at intermediate scales the energy cascades via the MHD non-linear coupling. The MHD cascade partially continues to small scales through the Hall term, but the small scales are dominated by the resistive dissipation; however, the pressure–strain effective dissipation takes place over a wide range of scales. The isotropised view of the KHM equation is semi-quantitatively similar to the 2D results (Hellinger et al. 2022). The oscillations of the compressive pressure–strain rate is accompanied by variations of the decay and the MHD cascade rates that are typically more important at larger scales.

The presence of the ambient magnetic field induces a strong anisotropy similar to that observed in 3D MHD simulations (Verdini et al. 2015; Montagud-Camps et al. 2022), even though the KHM equation is only very weakly dependent on B₀ (as clearly seen in its equivalent form, Eq. (3) of Hellinger et al. 2021a). The turbulent processes have similar anisotropic behaviour; they are weakly anisotropic on large scales (owing to the isotropic initialisation), and the anisotropy gets stronger towards smaller scales. The energy-containing range where the fluctuating energy decays is weakly anisotropic and importantly varies with time due to the compressive effects. The MHD cascade rate tends to be constant within an anisotropic range of scales. The Hall cascade rate is nearly stationary and strongly anisotropic; the onset of the Hall effect shifts to larger scales at quasi-parallel angles (with respect to the ambient magnetic field) compared to the quasi-perpendicular angles. The resistive dissipation is also nearly stationary and exhibits an anisotropy similar to that of the Hall cascade rate. The pressure–strain coupling is also anisotropic, works on a relatively wide range of scales, and has a weak anisotropy on large scales and strong one on small scales.

In the isotropised version of the KHM equation we have a relatively clear interpretation of the different terms. The cascade rates K_MHD and K_Hall represent the (magnetic+kinetic energy) cascade rates and ∂_tS, D, and Ψ have cumulative properties (cf. Hellinger et al. 2021a) as can be seen when comparing the KHM results to the results of the isotropic spectral transfer analysis (see Appendix A). The interpretation of the anisotropic terms is less clear. The behaviours of ∂_tS, D, and Ψ suggest some sort of cumulative properties. We tried a few spectral transfer analyses with some anisotropic filters (cylindrical or ellipsoid) and found only a qualitative similarity to the anisotropic KHM results. The cascade rates K_MHD and K_Hall seem to have the tendency to be constant in a region that is anisotropically distributed in the (l_⊥, l_∥). However, they are just scalars and indicate nothing about the direction of the cascading energy. The third-order structure function vectors in the usual form of the KHM equation possibly contain information about the direction of the cascade (cf. Verdini et al. 2015), but this assumption is not physically justified. Our simulation results, and previous results as well, indicate that turbulence develops in the inertial range K(l_⊥, l_∥) = Q with (generally anisotropic) outer scales connected with the initial conditions (or injection) properties. The turbulent cascade develops an important anisotropy that induces an anisotropy at the dissipation scales, and Hall scales as well.

The KHM equation is nicely satisfied in the simulation indicating that the kinetic plasma system and different turbulent processes are clearly numerically resolved. However, the simulation box size is not well separated from the ion characteristic scales owing to the computational constraints. Moreover, the ion pressure–strain coupling sets in at relatively large scales (∼10d_i) compared to the ion characteristic scales d_i and ρ_i. This is likely a reason for the fast onset and strong oscillations of the compressive pressure–strain rate, the decay, and the MHD cascade rate. We expect that for a larger system (and larger driving scales) the temporal variations would be reduced. This would justify our use of time-averaging to improve the statistics of turbulent fluctuations and to reduce the compressibility effects. More work is needed with a better separation between the injection (energy containing) range and the scale of the pressure–strain coupling. We presented one case study of plasma turbulence; more simulations are necessary to assess the role of the amplitude of turbulent fluctuations δB/B₀ and plasma betas, as well as the cross-helicity (Montagud-Camps et al. 2022).

Acknowledgments

This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254).

References

Appendix A: Isotropic spectral transfer

Another way to analyse the scale dependence of turbulence and its processes is spectral (Fourier) decomposition (cf. Mininni et al. 2007; Grete et al. 2017). To characterise turbulence it is possible to use isotropic low-pass filtered kinetic+magnetic energy (i.e., the energy in modes with wave-vector magnitudes smaller than or equal to k; cf. Hellinger et al. 2021b)

$$\begin{array}{cc}\hfill E_k& =\frac{1}{2}{\displaystyle \sum _{|{\mathit{k}}^{\prime }|\le k}(|\widehat{\mathit{w}}{|}^2+{|\widehat{\mathit{B}}|}^2),}\hfill \end{array}$$(A.1)

where the wide hat denotes the Fourier transform. For E_k we obtained the following dynamic equation

$$\begin{array}{c}\hfill {\partial }_t{E}_k+S_{\text{MHD}k}+S_{\text{Hall}k}=-{\mathrm{\Psi }}_k-D_k,\end{array}$$(A.2)

where S​_MHDk and S​_Hallk represent the MHD and Hall energy transfer rates, respectively; Ψ_k describes the pressure–strain effect; and D_k is the resistive dissipation rate for modes with wave-vector magnitudes smaller than or equal to k. They can be expressed as

$$\begin{array}{cc}\hfill S_{\text{MHD}k}& =\mathfrak{R}{\displaystyle \sum _{|{\mathit{k}}^{\prime }|\le k}[{\widehat{\mathit{w}}}^{\ast }·\widehat{(\mathit{u}·\mathbf{\nabla })\mathit{w}}+\frac{1}{2}{\widehat{\mathit{w}}}^{\ast }·\widehat{\mathit{w}(\mathbf{\nabla }·\mathit{u})}}\hfill \\ \hfill & -{\widehat{\mathit{w}}}^{\ast }·\widehat{{\rho }^{-1/2}\mathit{J}\times \mathit{B}}-{\widehat{\mathit{B}}}^{\ast }·\widehat{\mathbf{\nabla }\times (\mathit{u}\times \mathit{B})}],\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill S_{\text{Hall}k}& =\mathfrak{R}{\displaystyle \sum _{|{\mathit{k}}^{\prime }|\le k}{\widehat{\mathit{B}}}^{\ast }·\widehat{\mathbf{\nabla }\times (\mathit{j}\times \mathit{B})},}\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_k& =\mathfrak{R}{\displaystyle \sum _{|{\mathit{k}}^{\prime }|\le k}{\widehat{\mathit{w}}}^{\ast }·\widehat{{\rho }^{-1/2}\mathbf{\nabla }·\mathbf{P}},}\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill D_k& =\eta {\displaystyle \sum _{|{\mathit{k}}^{\prime }|\le k}|{\mathit{k}}^{\prime }{|}^2{|\widehat{\mathit{B}}|}^2.}\hfill \end{array}$$(A.6)

In these expression the asterisk denotes the complex conjugate, and ℜ denotes the real part. For fully developed turbulence S​_MHDk and S​_Hallk are the MHD and Hall cascade rates, respectively.

Similarly to the KHM case, we defined the validity test of the spectral transfer (ST) equation, Eq. (A.2), as

$$\begin{array}{c}\hfill O_k={\partial }_t{E}_k+S_{\text{MHD}k}+S_{\text{Hall}k}+{\mathrm{\Psi }}_k+D_k.\end{array}$$(A.7)

Figure A.1 shows the ST validity test O_k and the contributing terms as functions of k averaged over the time interval 200 ≤ tΩ_i ≤ 300.

Fig. A.1. Spectral transfer analysis, the validity test of the isotropic ST Ok, Eq. (A.7), as a function of k (averaged over the time interval 200 ≤ tΩi ≤ 300) (black line) along with the different contributing terms: decay rate ∂tEk (blue), MHD cascade rate S​MHDk (green), Hall cascade rate S​Hallk (orange), resistive dissipation rate Dk (red), and pressure–strain rate Ψk (magenta). All the quantities are given in units of the total (effective) heating rate Q.

Figure A.2 displays a direct comparison between the ST and KHM contributing terms as a function of k (through the inverse proportionality $kl=\sqrt{3}$. Figure A.2 shows that the MHD and Hall cascade rates in both the approaches are comparable

$$\begin{array}{c}\hfill S_{\text{MHD}k}\simeq K_{\mathrm{MHD}},S_{\text{Hall}k}\simeq K_{\mathrm{Hall}}.\end{array}$$(A.8)

Fig. A.2. Comparison between the ST and KHM approaches. The solid lines show the contributing terms of the ST equation as a function of k: (blue) the decay rate ∂tEk, (green) the MHD cascade rate SMHDk, (orange) the Hall cascade rate SHallk, (red) the dissipation rate Dk − Q, and (magenta) the pressure–strain rate Ψk − ψ. The dashed lines show the corresponding results of the KHM equation function of k through the relation $l=\sqrt{3}/k$: (blue) the decay rate −∂tS − Q, (green) the MHD cascade rate KMHD, (orange) the Hall cascade rate KHall, (red) the dissipation rate −D, and (magenta) the pressure–strain rate −Ψ. All the quantities are given in units of the total (effective) dissipation rate Q.

The MHD cascade term dominates at scales where the total power spectrum is close to a k^−5/3 power law. The other terms in the ST and KHM approaches are complementary:

$$\begin{array}{c}\hfill {\partial }_t[E_k+S(l)]\simeq -Q,D_k+D(l)\simeq Q_{\eta },{\mathrm{\Psi }}_k+\mathrm{\Psi }(l)\simeq \psi .\end{array}$$(A.9)

This means that ∂_tE_k represents the rate of change of the kinetic energy at scales with wave-vector magnitudes smaller than or equal to k, whereas ∂_tS represents approximately the rate for wave-vector magnitudes larger than k (with $l\simeq \sqrt{3}/k$). Similarly, D_k (Ψ_k) is the dissipation (pressure–strain) rate on scales ≤k, whereas D(l) (Ψ(l)) represents the corresponding rate on scales > k.

All Figures

Fig. 1. Evolution of different quantities as a function of time: (a) Relative changes in the kinetic energy ΔEkin (dashed line), magnetic energy ΔEmag (solid line), internal energy ΔEint (dash-dotted line), and total energy ΔEtot (dotted line); (b) resistive dissipation rate Qη (dashed line), incompressive pressure–strain effective dissipation rate ψi (solid line), and its compressive counterpart ψc (dotted line). These rates are normalised to the total effective dissipation rate Q = Qη + ψ (averaged over 200 ≤ tΩi ≤ 300).

In the text

	Fig. 2. Omnidirectional power spectral densities of the magnetic field B, PB (red), and the compensated proton velocity field w, Pw (blue), and their sum P (black) as a function of k normalised to di at $t=250{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$. The dotted line shows a spectrum ∝k−5/3 for comparison.
In the text

	Fig. 3. Schematic view of the theoretical expectations of a simplified system 0 = −∂tS + K − D. The decay term −∂tS (blue), the cascade rate K (green), and the dissipation rate −D (red) are given as a function of l and are normalised to the dissipation rate Q.
In the text

Fig. 4. Isotropised KHM analysis averaged results over a period of pressure–strain oscillations (200 ≤ tΩi ≤ 300): (a) Validity test O (black line) as a function of l along with combined dissipative terms −(D + Ψ) (red), combined cascade rates KMHD + KHall (green), and the decay rate −∂tS (blue) in a format similar to Fig. 3. (b) Different separate contributing terms to the KHM equation: the decay rate −∂tS (blue), the MHD cascade rate KMHD (green), the Hall cascade rate KHall (orange), the resistive dissipation rate −D (red), and the pressure–strain (effective dissipation) rate −Ψ (magenta). The dashed magenta curve displays the compressive component of −Ψ. All the quantities are normalised to the time-averaged effective total dissipation rate Q.

In the text

	Fig. 5. Evolution of isotropised KHM results. Shown are the different KHM terms as a function of time t and l: (a) decay rate ∂tS, (b) MHD cascade rate KMHD, (c) Hall cascade rate KHall, (d) resistive dissipation rate D, (e) compressive pressure–strain rate Ψc, and (f) incompressive pressure–strain rate Ψi. All the quantities are normalised to the effective total dissipation rate Q.
In the text

	Fig. 6. Color scale plots of gyrotropised power spectral densities of (a) the magnetic energy PB, and (b) the kinetic energy Pw as a function of k⊥ and k∥ (normalised to di).
In the text

Fig. 7. Anisotropic KHM results, labelled iso-contours of different quantities as a function of the separation scales l⊥ and l∥ (averaged over the time period 200 ≤ tΩi ≤ 300): (a) Second-order structure function S, (b) Decay rate ∂tS, (c) MHD cascade rate KMHD, (d) Hall cascade rate KHall, (e) pressure–strain (effective dissipation) term Ψ, and (f) resistive dissipation term D. The dotted line in panel a displays a circle $l_{\perp }^2+l_{\Vert }^2=\text{const}.$ as a measure of anisotropy. The dotted lines in panel c denote non-positive values. The KHM quantities are normalised to the (time-averaged) effective total dissipation rate Q.

In the text

Fig. A.1. Spectral transfer analysis, the validity test of the isotropic ST Ok, Eq. (A.7), as a function of k (averaged over the time interval 200 ≤ tΩi ≤ 300) (black line) along with the different contributing terms: decay rate ∂tEk (blue), MHD cascade rate S​MHDk (green), Hall cascade rate S​Hallk (orange), resistive dissipation rate Dk (red), and pressure–strain rate Ψk (magenta). All the quantities are given in units of the total (effective) heating rate Q.

In the text

Fig. A.2. Comparison between the ST and KHM approaches. The solid lines show the contributing terms of the ST equation as a function of k: (blue) the decay rate ∂tEk, (green) the MHD cascade rate SMHDk, (orange) the Hall cascade rate SHallk, (red) the dissipation rate Dk − Q, and (magenta) the pressure–strain rate Ψk − ψ. The dashed lines show the corresponding results of the KHM equation function of k through the relation $l=\sqrt{3}/k$: (blue) the decay rate −∂tS − Q, (green) the MHD cascade rate KMHD, (orange) the Hall cascade rate KHall, (red) the dissipation rate −D, and (magenta) the pressure–strain rate −Ψ. All the quantities are given in units of the total (effective) dissipation rate Q.

In the text