Ion-scale Transition of Plasma Turbulence: Pressure–Strain Effect

Turbulence in the weakly compressible solar wind plasma exhibits a clear transition at ion scales (Bruno & Carbone 2013). At large scales, magnetic power spectra of time series observed in situ have typically a power-law dependence on the frequency with a spectral index close to the Kolmogorov −5/3 phenomenological prediction for hydrodynamic turbulence. Around scales corresponding to ion characteristic scales (the ion gyroradius and inertial length) the power spectra steepen. This steepening was initially regarded as a signature of the dissipation onset (Leamon et al. 1998). However, the Hall term starts to play on similar scales and leads also to a spectral steepening without necessarily implying the presence of some energy dissipation (Ghosh et al. 1996; Galtier 2006; Papini et al. 2019).

One way to discern and quantify different turbulent processes is the Kármán–Howarth–Monin (KHM) equation (de Kármán & Howarth 1938; Monin & Yaglom 1975; Frisch 1995) that connects the energy decay/injection with its cascade and dissipation. Recently, the incompressible version of the KHM equation for the Hall–magnetohydrodynamic (MHD) approximation (Politano & Pouquet 1998; Galtier 2008; Hellinger et al. 2018; Ferrand et al. 2019) was used to study the ion transition in simulations as well as observations (Hellinger et al. 2018; Bandyopadhyay et al. 2020; Adhikari et al. 2021). These results indicate that at ion scales there is indeed a transition from an MHD- to a Hall-dominated turbulent cascade. However, through this transition the total cascade rate is observed to decrease and this suggests some sort of dissipation. Consequently, based on the previous works above, the ion transition is likely a combination of the Hall physics onset and dissipation.

Theoretical analyses (Yang et al. 2017, and references therein) show that one possible channel that leads to energy exchanges between the magnetic + particle kinetic energy and the particle internal energies in collisionless plasmas is the pressure–strain effect, a generalization of the pressure–dilation effect. Results of kinetic simulations (Yang et al. 2019; Matthaeus et al. 2020) indicate that this effect can act as an effective dissipation and may explain the decrease of the energy cascade rate observed at ion scales. Moreover, numerical simulations and in situ observations suggest that the pressure–strain channel is likely responsible for the correlations between particle-velocity-field gradients and temperatures (Franci et al. 2016a; Parashar & Matthaeus 2016; Yang et al. 2019; Pezzi et al. 2021; Yordanova et al. 2021).

In this paper we revisit the work of Hellinger et al. (2018) and analyze the kinetic simulation results using both the incompressible and compressible versions of the KHM equation in order to test the validity of these approximations. We use the compressible KHM equation derived by Hellinger et al. (2021a) because, motivated by the previous works of Yang et al. (2017) and Matthaeus et al. (2020), we want to determine effects of the pressure–strain coupling while the alternative approaches (Andrés et al. 2018, and references therein) assume a scalar pressure along with some particular closure, thus preventing their extension to the weakly collisional case with a tensor description of the particle pressure.

This paper is organized as follows: In Section 2 we present spectral properties of three two-dimensional (2D) hybrid simulations. In Section 3 we analyze the hybrid simulations using the incompressible KHM equation. In Section 4 we extend the analysis to the compressible KHM equation. In Section 5 we summarize and discuss the presented results.

Here we analyze 2D hybrid simulations of decaying plasma turbulence using the KHM equations. In the hybrid approximation, ions are described by a particle-in-cell model whereas electrons are a massless, charge-neutralizing fluid (Matthews 1994). We use the 2D version of the code Camelia ⁸ using a simulation setup that is similar to that of Franci et al. (2015).

We investigate properties of three hybrid simulations with parameters similar (but not identical) to those in Hellinger et al. (2018); in contrast to this work, we try here to use the same/similar parameters for all the runs if possible, see Table 1. Protons are initially isotropic with different values of β_i (run 1: 0.1, run 2: 0.5, and run 3: 2.5). We consider a 2D domain (x,y) of size 2048 × 2048 grid points and resolution Δx = Δy = d_i /16 (for runs 1 and 2) and d_i /8 (for run 3). Here d_i denotes the ion inertial length and β_i stands for the ion beta, i.e., the ratio between the ion and 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 and perpendicular to the simulation domain is present whereas neutralizing electrons are assumed to be isotropic and isothermal. Furthermore, we set the electron beta (i.e., the ratio between electron and magnetic pressures) as equal to the initial value of the ion beta (β_e = β_i ). The system is perturbed with an isotropic 2D spectrum of modes with random phases, linear Alfvén polarization (δ B ⊥ B ₀), and vanishing correlation between magnetic field and velocity fluctuations. These modes are in the range $k\leqslant 0.22{d}_{i}^{-1}$ (for runs 1 and 2) and $k\leqslant 0.2{d}_{i}^{-1}$ (for run 3) and have a flat one-dimensional power spectrum with rms fluctuations δ B = 0.25. The time step is ${\rm{\Delta }}t=0.005{{\rm{\Omega }}}_{i}^{-1}$ (for runs 1 and 2) and $0.01{{\rm{\Omega }}}_{i}^{-1}$ (for run 3) 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 = 8192, and a small resistivity $\eta =4\times {10}^{-4}{\mu }_{0}{v}_{A}^{2}/{{\rm{\Omega }}}_{i}$ is used to avoid energy accumulation at the smallest scales (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 beyond the time when it becomes quasi-stationary (Mininni & Pouquet 2009; Servidio et al. 2015); henceforth, we analyze properties of plasma turbulence at such times ${t}_{d}=286{{\rm{\Omega }}}_{i}^{-1}$, $290{{\rm{\Omega }}}_{i}^{-1}$, and $299{{\rm{\Omega }}}_{i}^{-1}$ for runs 1, 2, and 3, respectively.

Figure 1 shows the power spectral density of the magnetic field for the three simulations. The simulated spectra exhibit two power laws with a smooth transition at ion scales (the so-called ion spectral break), whose shape and position depend on the plasma beta: its scale is close to d_i for small betas whereas in high-beta plasmas it is around ρ_i (Franci et al. 2016b). Here ρ_i denotes the ion gyroradius. The magnetic power spectral slopes on large scales are about −5/3, in the subion range the spectrum steepens, with slopes about −3.5, −3, and −2.9 in the three simulations (see Franci et al. 2016b). On the other hand, the magnetic compressibility increases at ion scales depending on the plasma beta (Matteini et al. 2020). The proton velocity fluctuating field has a power spectrum on large scales with a similar slope and decouples from the magnetic fluctuations around the proton gyroscales. The subion velocity fluctuations have a limited scale range before reaching the noise level so that it is difficult to distinguish between an exponential and a power-law dependence in the subion range.

Zoom In Zoom Out Reset image size

We start with the incompressible, constant-density, inviscid (there is no explicit viscosity in the hybrid model), but resistive Hall–MHD approximation. In this case we have the energy (per unit mass) budget equation for the kinetic+magnetic energy

$$\begin{eqnarray}&&\frac{\partial }{\partial t}\lt\frac{1}{2}\left(| {\boldsymbol{u}}{| }^{2}+| {\boldsymbol{b}}{| }^{2}\right)\gt-\epsilon ,\end{eqnarray} \tag{ 1 }$$

where u is the velocity field, b is the magnetic field in the Alfvén units ( b = B /ρ^1/2, B being the magnetic field and ρ the plasma density assumed to be constant; here we assume SI units except for the magnetic permeability μ₀, which is set to one), and

$$\begin{eqnarray}&&\epsilon =\eta \langle {\boldsymbol{\nabla }}{\boldsymbol{b}}:{\boldsymbol{\nabla }}{\boldsymbol{b}}\rangle \end{eqnarray} \tag{ 2 }$$

is the resistive dissipation rate (per unit mass). Following Hellinger et al. (2018) we define the effective dissipation/cascade rate (per unit mass) ${\epsilon }^{* }$ given by the corresponding KHM equation (the superscript (i) will henceforth denote incompressible, constant-density quantities):

$$\begin{eqnarray}&&{\epsilon }^{* }=-\frac{1}{4}\frac{\partial {S}^{(i)}}{\partial t}+{K}_{\mathrm{MHD}}^{(i)}+{K}_{{\rm{H}}}^{(i)}+\frac{\eta }{2}{\rm{\Delta }}{S}_{b}^{(i)},\end{eqnarray} \tag{ 3 }$$

where ${S}_{b}^{(i)}$, ${S}_{u}^{(i)}$, and ${S}^{(i)}={S}_{b}^{(i)}+{S}_{u}^{(i)}$ are the second-order structure functions:

$$\begin{eqnarray}&&{S}_{b}^{(i)}=\langle {\left|\delta {\boldsymbol{b}}\right|}^{2}\rangle ,{S}_{u}^{(i)}=\langle {\left|\delta {\boldsymbol{u}}\right|}^{2}\rangle ,\end{eqnarray} \tag{ 4 }$$

${K}_{\mathrm{MHD}}^{(i)}$, and ${K}_{{\rm{H}}}^{(i)}$ are given by

$$\begin{eqnarray}&&{K}_{\mathrm{MHD}}^{(i)}=-\displaystyle \frac{1}{4}{\boldsymbol{\nabla }}\cdot {{\boldsymbol{Y}}}^{(i)},\ \mathrm{and}\ {K}_{{\rm{H}}}^{(i)}=-\displaystyle \frac{1}{4}{\boldsymbol{\nabla }}\cdot {{\boldsymbol{H}}}^{(i)},\end{eqnarray} \tag{ 5 }$$

respectively, involving the third structure functions:

$$\begin{eqnarray}&&{{\boldsymbol{Y}}}^{(i)}=\left\langle \delta {\boldsymbol{u}}\left({\left|\delta {\boldsymbol{u}}\right|}^{2}+{\left|\delta {\boldsymbol{b}}\right|}^{2}\right)-2\delta {\boldsymbol{b}}\left(\delta {\boldsymbol{u}}\cdot \delta {\boldsymbol{b}}\right)\right\rangle ,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\boldsymbol{H}}}^{(i)}=\left\langle \delta {\boldsymbol{b}}\left(\delta {\boldsymbol{b}}\cdot \delta {\boldsymbol{j}}\right)-\displaystyle \frac{1}{2}\delta {\boldsymbol{j}}{\left|\delta {\boldsymbol{b}}\right|}^{2}\right\rangle .\end{eqnarray} \tag{ 7 }$$

Here δ denotes the increment for the spatial separation l , δ u = u ( x + l ) − u ( x ), etc., j is the electric current in the Alfvén units ( j = J /(en), J being the electric current, e the proton charge, and n the proton number density). In Equation (3), ${S}_{u}^{(i)}$ and ${S}_{b}^{(i)}$ represent the separation-scale distribution of the kinetic and magnetic energy (per unit mass), respectively, ${K}_{\mathrm{MHD}}^{(i)}$ and ${K}_{{\rm{H}}}^{(i)}$ are the MHD and Hall cascade rates, respectively, and $\eta {\rm{\Delta }}{S}_{b}^{(i)}/2$ describes the resistive dissipation.

Equation (3) constitutes an energy-per-mass budget relationship. It describes the balance between injection, energy cascade flux, and dissipation; the KHM equation corresponds to the case when ${\epsilon }^{* }=\epsilon =\mathrm{const}$. We can now directly verify the predicted constancy of the right-hand side of Equation (3) in the simulations, and compare the relative importance of the different terms across the separation scales. Figure 2 shows the results of such analysis, displaying the effective cascade rate ${\epsilon }^{* }$ and its contributing terms as functions of the scale separation l for the three simulations. The different terms are calculated at the times t_d ; ∂S/∂t is estimated by the finite difference (S(t_d + Δt) − S(t_d ))/Δt using the simulation (particle) time step.

Zoom In Zoom Out Reset image size

Figure 2 demonstrates that the effective cascade rate in the three simulations ${\epsilon }^{* }$ is constant on large scales but larger than the expected resistive dissipation rate . The cascade rate ${\epsilon }^{* }$ decreases at ion scales and this transition shifts to larger scales for larger β_i (and larger ρ_i ). This behavior was already seen in Hellinger et al. (2018) but in this work the decrease was underestimated since the Hall terms had twice the correct value (see Ferrand et al. 2019). These properties indicate the existence of an additional dissipation channel not captured in this analysis.

We now relax the incompressibility and constant-density assumption (but we still neglect the viscosity) and, moreover, we consider a weakly collisional plasma where the pressure is described as a pressure tensor ( P ). The kinetic+magnetic energy budget equation reads then

$$\begin{eqnarray}&&\frac{\partial }{\partial t}\lt\frac{1}{2}\left(| {\boldsymbol{w}}{| }^{2}+| {\boldsymbol{B}}{| }^{2}\right)\gt=-Q,\end{eqnarray} \tag{ 8 }$$

where we include the (square root of the) plasma density into the velocity field ${\boldsymbol{w}}=\sqrt{\rho }{\boldsymbol{u}}$ (Kida & Orszag 1990), and where we define the total effective dissipation rate (per unit volume here in contrast to the incompressible case)

$$\begin{eqnarray}&&Q={Q}_{\eta }+{Q}_{{\rm{\Sigma }}}\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}&&{Q}_{\eta }=\eta \langle {\boldsymbol{\nabla }}{\boldsymbol{B}}:{\boldsymbol{\nabla }}{\boldsymbol{B}}\rangle =\eta \langle {J}^{2}\rangle \end{eqnarray} \tag{ 10 }$$

being the resistive dissipation rate (simply related to the resistive dissipation rate per unit mass as Q_η = 〈ρ〉; here 〈ρ〉 is the background/average plasma density needed to relate the normalizations per unit mass and per unit volume) and

$$\begin{eqnarray}&&{Q}_{{\rm{\Sigma }}}=-\langle {\boldsymbol{P}}:{\boldsymbol{\Sigma }}\rangle \end{eqnarray} \tag{ 11 }$$

representing the pressure–strain coupling (Yang et al. 2017; Σ being the velocity–strain tensor, Σ = ∇ u ). To generalize Equation (3), we consider the compressible KHM equation (Equation (4) of Hellinger et al. 2021a); we express the resistive dissipation rate as a function of the remaining terms and add the pressure–strain term Q_Σ to both sides and we define an effective cascade rate Q* as

$$\begin{eqnarray}&&{Q}^{* }=-\displaystyle \frac{1}{4}\displaystyle \frac{\partial S}{\partial t}+{K}_{\mathrm{MHD}}+{K}_{{\rm{H}}}+{\rm{\Psi }}+\displaystyle \frac{\eta }{2}{\rm{\Delta }}{S}_{B}.\end{eqnarray} \tag{ 12 }$$

Here the second-order structure functions are given as

$$\begin{eqnarray*}&&{S}_{w}=\left\langle | \delta {\boldsymbol{w}}{| }^{2}\right\rangle ,\ {S}_{B}=\left\langle | \delta {\boldsymbol{B}}{| }^{2}\right\rangle ,\ \mathrm{and}\ S={S}_{w}+{S}_{B},\end{eqnarray*}$$

and the following quantities:

$$\begin{eqnarray}\begin{array}{rcl}{K}_{\mathrm{MHD}} & = & -\frac{1}{4}{\boldsymbol{\nabla }}\cdot {\boldsymbol{Y}}-\frac{1}{4}R+\frac{1}{2}{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{B}}\times {\boldsymbol{J}}\right],\\ {K}_{{\rm{H}}} & = & -\frac{1}{4}{\boldsymbol{\nabla }}\cdot {\boldsymbol{H}}+\frac{1}{4}{C}_{\rho }\left[{\boldsymbol{B}}\times {\boldsymbol{j}},{\boldsymbol{J}}\right]\end{array}\end{eqnarray} \tag{ 13 }$$

are connected, similarly to the previous case, with the third-order structure functions:

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{Y}} & = & \left\langle \delta {\boldsymbol{u}}\left(| \delta {\boldsymbol{w}}{| }^{2}+| \delta {\boldsymbol{B}}{| }^{2}\right)-2\delta {\boldsymbol{B}}\left(\delta {\boldsymbol{u}}\cdot \delta {\boldsymbol{B}}\right)\right\rangle ,\\ {\boldsymbol{H}} & = & \left\langle \delta {\boldsymbol{B}}\left(\delta {\boldsymbol{j}}\cdot \delta {\boldsymbol{B}}\right)-\displaystyle \frac{1}{2}\delta {\boldsymbol{j}}| \delta {\boldsymbol{B}}{| }^{2}\right\rangle .\end{array}\end{eqnarray} \tag{ 14 }$$

The third-order structure functions Y and H represent the compressible generalizations of Y ⁽ⁱ⁾ and H ⁽ⁱ⁾ (note that H = 〈ρ〉 H ⁽ⁱ⁾), whereas R is a compressible term that does not seem to be easily expressible as a structure function but can be given as

$$\begin{eqnarray}&&R=\lt\delta {\boldsymbol{w}}\cdot \left({\theta }^{^{\prime} }{\boldsymbol{w}}-\theta {{\boldsymbol{w}}}^{{\prime} }\right)\gt,\end{eqnarray} \tag{ 15 }$$

where θ is the dilatation field,

$$\begin{eqnarray}&&\theta ={\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 16 }$$

and where the primes denote quantities evaluated at ${{\boldsymbol{x}}}^{{\prime} }\,={\boldsymbol{x}}+{\boldsymbol{l}}$. Equation (13) defining K_MHD and K_H also involves "correction" terms related to the density variations given by

$$\begin{eqnarray*}&&{C}_{\rho }\left[{\boldsymbol{a}},{\boldsymbol{b}}\right]=\left\langle \left(\displaystyle \frac{{\rho }^{{\prime} }}{\rho }-1\right){{\boldsymbol{a}}}^{{\prime} }\cdot {\boldsymbol{b}}+\left(\displaystyle \frac{\rho }{{\rho }^{{\prime} }}-1\right){\boldsymbol{a}}\cdot {{\boldsymbol{b}}}^{{\prime} }\right\rangle .\end{eqnarray*}$$

Similarly to the incompressible Equation (3), in Equation (12) S_w and ${S}_{B}=\langle \rho \rangle {S}_{b}^{(i)}$ represent the separation-scale distribution of the kinetic and magnetic energy, K_MHD and K_H are the MHD and Hall cascade rates, respectively, and ηΔS_B /2 describes the resistive dissipation. Furthermore, Ψ represents the pressure–strain effect,

$$\begin{eqnarray}&&{\rm{\Psi }}=-\displaystyle \frac{1}{2}\left\langle \delta {\boldsymbol{w}}\cdot \delta \left(\displaystyle \frac{{\boldsymbol{\nabla }}\cdot {\boldsymbol{P}}}{\sqrt{\rho }}\right)\right\rangle +{Q}_{{\rm{\Sigma }}}.\end{eqnarray} \tag{ 17 }$$

This term is different from that defined in Hellinger et al. (2021a), here we have added the pressure–strain rate Q_Σ.

Equation (12) describes an energy budget relationship and the KHM equation corresponds to the equality ${Q}^{* }=Q=\mathrm{const}$. Figure 3 shows the test of the KHM equation, the cascade rate Q*, and its contributing terms as functions of the separation scale l in the three simulations. Figure 3 demonstrates that the compressibility and inclusion of the pressure–strain effect significantly improves the conservation of the effective cascade/dissipation rate. Q* is relatively constant in the three simulations, and the decrease of the MHD + Hall cascade rate observed at the transition from the MHD to subion scales is compensated by the resistive and pressure–strain coupling. The contribution from this pressure–strain coupling term is of similar amplitude as the resistive one at low and moderate beta, while it dominates in the high-beta case of run 3. Besides the constancy, Q* is close to its expected value Q; the relative error is small ∣Q^*∣/Q ≲ 0.1 and is likely related to numerical limitation of the code.

Zoom In Zoom Out Reset image size

Comparing more in detail individual contributions in the two approaches, we see that the compressibility and density variations are not important in the current simulations. The Hall contributions are almost the same, $\langle \rho \rangle {K}_{{\rm{H}}}^{(i)}\simeq {K}_{{\rm{H}}}$, and the MHD contributions are close to each other, $\langle \rho \rangle {K}_{\mathrm{MHD}}^{(i)}\sim {K}_{\mathrm{MHD}}$; the relative difference is about 10% for β_i = 0.1 (mostly owing to the compressible term R) and decreases for higher betas. The main difference between the presentations of these terms in Figures 2 and 3 is owing to the different normalizations ( versus Q); this normalization has actually a highly physical meaning that expresses the inertial-range picture relating the cascade rate to the dissipation one. In order to have the cascade rate equal the dissipation one, one has to include the pressure–strain contribution to the dissipation rate. We conclude that the substantial difference between Equations (3) and (12) for the three simulations investigated in this paper is owing to the inclusion of the pressure–strain effect.

In this paper we investigated the transition from MHD to subion scales in 2D hybrid simulations using the KHM equations. We analyzed three 2D hybrid simulations with parameters similar (but not identical) to those in Hellinger et al. (2018). We observe qualitatively the same results as in Hellinger et al. (2018), the incompressible, constant-density KHM equation shows that the effective cascade/dissipation rate decreases at subion scales. Note that the incompressible Hall cascade rate in Hellinger et al. (2018) is twice the correct value (see Ferrand et al. 2019) so that the decrease of the effective incompressible cascade rate at subion scales is underestimated there. The decrease of the incompressible rate at subion scales is in agreement with other previous simulation and observation works (Bandyopadhyay et al. 2020; Adhikari et al. 2021).

The compressible version of the KHM equation (Hellinger et al. 2021a) that includes the pressure–strain coupling (Yang et al. 2017) exhibits a good conservation property and the effective cascade/dissipation rate is constant; the decrease of the cascade rate at ion scales is compensated by the resistive effects, and, more importantly, the pressure–strain coupling. The effective cascade/dissipation rate is close to the expected value. We observe a small but nonnegligible discrepancy that is likely connected with numerical issues of the particle-in-cell, finite difference scheme. The discrepancy between the prediction of the KHM equation and the simulation results depends on many numerical as well as physical parameters. In fact, the KHM equation may serve as a test of numerical codes. We were not able to discern and analyze all the relevant parameters. Additional simulations show that the error decreases with an increasing number of particles per cell and decreasing beta. This indicates that the numerical noise due to the limited number of particles per cell contributes to the error; the finite difference scheme for the resistivity is a possible source of additional errors, and the pressure–strain coupling is prone to be affected by the numerical noise.

In a more general context, our study suggests that, for quite a wide range of plasma and turbulence parameters, the actual compressibility and density variations are weak; one can use an incompressible approach but pressure–strain effect needs to be retained. The same probably applies for the weakly compressible solar wind (Zank et al. 2017). The spatial scale of the pressure–strain onset increases with β_i (and ρ_i ). Our current results do not, however, show a clear connection between this onset scale and ρ_i and/or d_i . The characteristic scales of the pressure–strain channel, P : Σ, are apparently difficult to determine (Del Sarto & Pegoraro 2018).

The pressure–strain effect is in principle reversible but, in contrast to fluid cases (Hellinger et al. 2021a, 2021b), it appears to work only in one direction and acts as an effective dissipation rate in the kinetic case. Our kinetic/hybrid KHM results are similar to the spatial filtering analysis of fully kinetic simulations (Matthaeus et al. 2020) and constitute a complementary and independent confirmation that pressure–strain coupling plays the role of a dissipation channel that appears at ion scales. More kinetic simulations are needed to discern its characteristic scales and the roles of different ion and electron species and their parameters; the KHM equation needs to be extended to the multifluid framework (Andrés et al. 2016). We also note that multiple different processes may contribute to the pressure–strain effect, such as quasilinear damping, magnetic reconnection, etc., and one has to go beyond the KHM equation to distinguish them.

We plan to continue this work using the spectral transfer analysis that allows a clearer connection between the different energy-transfer channels and the power spectra (see Papini et al. 2021). We also plan to extend this analysis to a three-dimensional geometry (see Verdini et al. 2015; Franci et al. 2018) to investigate the anisotropy of the energy cascade and dissipation including the pressure–strain effect. Results obtained in this work by means of numerical simulations need to be complemented by in situ observations. In this respect, it is already challenging to extend the incompressible KHM equation to the Hall regime (Bandyopadhyay et al. 2020). However, multipoint spacecraft observations (by, e.g., the Magnetospheric Multiscale Mission) may possibly be used to measure the compressible Hall–MHD in situ, also including the pressure–strain effect (Bandyopadhyay et al. 2021).

This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility. The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grant Nos. ST/P002307/1 and ST/R002452/1 and STFC operations grant No. ST/R00689X/1. This work also used computing resources provided by STFC DiRAC HPC Facility at Durham (grant Nos. ST/P002293/1, ST/R002371/1, ST/S002502/1, ST/R000832/1) for project "dp170" and by Cineca and INAF (Accordo Quadro MoU Nuove frontiere in Astrofisica) for project "INA20_C6A55". L.F. is supported by the STFC grant No. ST/T00018X/1.