Spacetime Hall-MHD Turbulence at Sub-ion Scales: Structures or Waves?

The physical mechanisms underlying the turbulent cascade in the solar wind and other space plasma environments below the proton characteristic scales remain largely unknown. In particular, it is hotly debated whether the energy cascade at sub-ion scales results from nonlinearly interacting wavelike fluctuations described in terms of low-frequency kinetic Alfvén waves (KAW; Sahraoui et al. 2009; Schekochihin et al. 2009; Salem et al. 2012; Chen et al. 2013b) and/or whistler waves, eventually complemented by wave–particle interactions (e.g., Landau damping, Sahraoui et al. 2010; Sulem et al. 2016), or from spatially localized highly intermittent coherent structures (e.g., current sheets and reconnection sites) where dissipation, heating, and cross-scale energy transfer are enhanced (Perri et al. 2012; Wan et al. 2012; Osman et al. 2014; Cerri & Califano 2017; Yang et al. 2017; Camporeale et al. 2018; Papini et al. 2020). It is difficult to address this question, since the intrinsic spatiotemporal multiscale nature of turbulence makes it difficult to identify and isolate the relevant mechanisms in place. At MHD scales and for moderate values of the plasma β, the turbulent eddies carry a significant part of the energy (Andrés et al. 2017). At sub-ion scales, early studies highlighting the spacetime Fourier structure of turbulence hint that the plasma dynamics is dominated by almost static low-frequency features, with no clear signature of waves (Parashar et al. 2010). Recently, several studies that exploit advanced techniques for the analysis of the spatial properties of turbulence provided evidence that the cross-scale energy transfer at sub-ion scales occurs preferentially in coherent structures through enhanced dissipation (Yang et al. 2017; Camporeale et al. 2018; Papini et al. 2020). On the one hand, such a picture is consistent with the high levels of intermittency measured in the solar wind (Osman et al. 2012; Wu et al. 2013). On the other hand, solar wind turbulence at sub-ion scales shown to have properties similar to those of KAWs, such as well-defined ratios between the power spectra of magnetic and density fluctuations (Salem et al. 2012; Chen et al. 2013b) and a reduced magnetic helicity (Matthaeus et al. 1982) of positive sign at high angles of propagation with respect to the mean magnetic field (Telloni & Bruno 2016; Woodham et al.2018).

In this Letter, we report results from a high-resolution spatiotemporal multiscale study of fully developed turbulence in a 2D Hall-magnetohydrodynamic (MHD) simulation, with the aim to assess whether nonlinear wavelike activity at sub-ion scales is energetically relevant to the plasma dynamics. Previous studies, comparing Hall-MHD models with hybrid (ion-kinetic and electron-fluid) models (Papini et al. 2019) and fully kinetic models (González et al. 2019), show that Hall-MHD is able to reproduce the main magnetic properties of plasma turbulence at sub-ion scales. In turn, hybrid kinetic models well reproduce in situ observations (Franci et al. 2020, 2020). The use of a 2D setup allows for a much higher resolution, compared to a full 3D one, while preserving many spectral and statistical properties of turbulence (Franci et al. 2018). It is true, however, that employing a 2D geometry prevents the onset of 3D intrinsic features of plasma turbulence (e.g., critically balanced turbulence; Goldreich & Sridhar 1995). The novelty of this work relies in the concurrent use of (i) a large high-resolution spatiotemporal data set and, (ii) Fast Iterative Filtering (FIF), a novel adaptive technique designed for the decomposition of nonlinear nonstationary signals (Cicone 2020; see Appendix A), which allows us to clearly isolate features in the time-frequency domain.

We performed a simulation of freely decaying Alfvénic turbulence on a box of size L_x × L_y = 128d_i × 128d_i , with 1024² gridpoints and a resolution Δx = Δy = d_i /8 (where d_i is the ion inertial length), by means of a fully compressible viscous-resistive nonlinear Hall-MHD pseudospectral code (see Appendix C and Papini et al. 2019). The system is initialized with a constant out-of-plane magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$ ($\hat{{\boldsymbol{z}}}$ defines the parallel direction) and a uniform plasma β = 2. Moreover, the xy-plane is filled with incompressible velocity ( u ) and magnetic ( B ) fluctuations in energy equipartition and perpendicular to B ₀. These fluctuations are random-phase sinusoidals of constant amplitude and with wavenumbers spanning − 8π/L_x(y) ≤ k_x(y) ≤ 8π/L_x(y). Their initial rms amplitude is u_rms = B_rms ≃ 0.24.

Our analysis begins at the maximum of turbulent activity at ${t}_{0}=195{{\rm{\Omega }}}_{i}^{-1}$ (where Ω_i is the ion-cyclotron (IC) frequency), corresponding to the peak of the rms of the current density (Mininni & Pouquet 2009). The spacetime analysis is performed in the time interval $[{t}_{0},{t}_{1}]=[195{{\rm{\Omega }}}_{i}^{-1},215{{\rm{\Omega }}}_{i}^{-1}]$, in which all fields are sampled with a cadence of $0.01{{\rm{\Omega }}}_{i}^{-1}$. The resulting fields are defined on a spatiotemporal grid of 1024 × 1024 × 2001 points. We separate the contribution of fast and slowly evolving features by means of a new approach that uses Fourier transform (FT) and FIF (Lin et al. 2009; Cicone et al. 2016) for the spatial and the temporal decomposition, respectively. We name such a mixed approach as the FTFIF decomposition. The novelty of our analysis is the use of FIF for the temporal decomposition. FIF is designed to decompose a nonstationary nonlinear signal into a set of intrinsic mode components (IMCs) oscillating around zero but with varying amplitude and frequency, plus a residual or trend. Such decomposition is adaptive, based on the local characteristic frequencies of the signal, and does not make any assumption on the shape of the signal to be extracted. As such, FIF offers many advantages over more traditional methods (e.g., FT and Wavelets) that require assumptions on the stationarity and/or linearity of the signal. Further details on the FIF technique and on our FTFIF analysis are found in Appendices A and B.

Figure 1 shows the k ω-diagram of the power spectrum (or power spectral density) P(k_⊥, ω) for both the parallel ( B _z ) and perpendicular ( B _⊥) magnetic field and for the perpendicular velocity fluctuations ( u _⊥), as obtained by the FTFIF decomposition. The k ω-diagram of a field is computed by first performing an FT in space. Then, a FIF decomposition is performed in time to obtain a set of IMCs for each Fourier mode. Finally, the power and average temporal frequency of each IMC is interpolated to a k ω-grid and summed over all the IMCs (see Appendix B). On average, more than 2,000,000 IMCs were extracted per field component.

Zoom In Zoom Out Reset image size

The k ω-diagram of u _⊥ (Figure 1(a)) shows clear signatures of wave activity. Two ridges emanating from the origin nicely follow the theoretical dispersion relations (Pucci et al. 2016) for fast/whistler (FW) waves and Alfvén/kinetic Alfvén waves (A/KAW) respectively. The dispersion relation for A/KAW is calculated by using the mean (79.3°) of the angle θ_Bk between the magnetic field and the xy-plane (where the k vectors lie). Such a ridge is broadened because of the excursion of the values of θ_Bk and β in the simulation at the maximum of the turbulent activity, which induces changes in the frequencies of the waves throughout the simulation domain. For FW, the best matching relation is obtained by using the mode (84.9°) of θ_Bk . This, together with the fact that the A/KAW ridge shows a wider broadening, suggests that FW are preferentially excited in regions with very small perpendicular magnetic fluctuations where the magnetic field is almost parallel to the mean field B ₀, while KAWs are more widely distributed. According to its polarizations, the FW ridge is visible also in the k ω-spectrum of B _z and ρ (not shown). The A/KAW ridge is not as visible, since it is partially swamped by the turbulent magnetic activity. In the same figure, we report the frequency ω_nl = 2π/τ_nl corresponding to the characteristic nonlinear turbulent time ${\tau }_{{nl}}({k}_{\perp })=1/({k}_{\perp }{\widetilde{u}}_{e}({k}_{\perp }))$, which estimates the rate at which turbulence evolves at a given scale ℓ = 2π/k_⊥ (Papini et al. 2019) (${\widetilde{u}}_{e}({k}_{\perp })$ being the amplitude of the electron fluid velocity at that scale). It is reasonable to assume that the area in the k ω-diagram below ω_nl is dominated by the turbulent dynamics whereas, above ω > ω_nl , features (such as the FW ridge and part of the A/KAW ridge) can evolve freely and independently from the main turbulent cascade.

To evaluate the energy contribution of FW and A/KAW, we calculated their k-power spectra by integrating in frequency P(k_⊥, ω) over the FW and the A/KAW ridge, respectively. The integration area of the FW ridge encompasses all the frequencies ω within 10% of the eigenfrequencies of the FW dispersion relation. Instead, the integration area of the A/KAW ridge was obtained by considering the excursion of the values of θ_Bk and β in the simulation domain, which affects the frequencies of KAWs throughout the domain, thus causing the broadening of the ridge. Results are shown in Figure 2 for u _⊥, B _⊥, and B _z (but similar results are obtained for the other fields). For all the fields, both the FW (in violet) and KAW (in orange) spectra are energetically negligible at sub-ion scales, being more than one order of magnitude smaller than the total power (in black). The A/KAW spectrum is not negligible at fluid-MHD scales, as expected in the case of Alfvénic turbulence.

Zoom In Zoom Out Reset image size

Although energetically irrelevant at sub-ion scales, both the FW and the A/KAW spectra appear to follow well-defined power laws. At around the spectral break (k_⊥ d_i ≃ 2.8, vertical dotted line) both the FW spectrum (in B _z ) and the A/KAW spectrum (in B _⊥ and B _z ) sharply transition to a well-defined −11/3 slope, which persists until dissipation scales are reached (at k_⊥ d_i ≃ 12, vertical dotted–dashed line). The FW spectrum of u _⊥ also shows a −11/3 power law, because of its polarization properties. The existence of a −11/3 power law at sub-ion scales is predicted by few Hall-MHD theoretical models (Krishan & Mahajan 2004; Galtier & Buchlin 2007; Meyrand & Galtier 2012; Schekochihin et al. 2019). Some of their underlying hypotheses, however, are inconsistent with what was observed here (e.g., the requirement of P _u (k_⊥) ≫ P _B (k_⊥) in Krishan & Mahajan 2004, Galtier & Buchlin 2007, and Schekochihin et al. 2019 is inconsistent with both the FW and A/KAW spectra). This, together with the fact that the whole FW ridge and the A/KAW ridge (although partially) are located at frequencies higher than the nonlinear-time frequency ω_nl (see Figure 1), possibly suggests that FW and KAWs may undergo a separate cascade independently from the main cascade. We remark that both the FW and the A/KAW spectra remain energetically irrelevant even if the integration area over the corresponding ridge is increased, although the −11/3 power law then tends to disappear.

Instead of being related to FW and/or KAW activity, the energy at sub-ion scales is concentrated at low frequencies. For instance, Figure 3 shows that the power spectrum of B _⊥ (solid black line) mainly results from the energetic contribution of slowly evolving structures and/or perturbations with temporal frequencies smaller than the IC frequency Ω_i (solid blue line). Going toward smaller scales (k_⊥ d_i ≳ 6), the contribution of medium frequencies (ω ≤ 3Ω_i , solid red curve) becomes important. The energy contained at higher frequencies (dashed red curve) is negligible.

Zoom In Zoom Out Reset image size

The simulation data show that the turbulent energy at kinetic scales is concentrated at low frequencies and it is not clearly related to KAW activity. At first glance, these results seem to disagree with previous studies (Chen et al. 2013a) which used the ratio between magnetic and density power spectra (the so-called KAW ratios) to detect the presence of KAW activity in the solar wind. We deem here that it may not be correct to use such ratios as a proxy for KAWs (as also noted by, e.g., Grošelj et al. 2019). To support such statement, in Figure 4(a) we report the ${{\rm{P}}}_{{{\boldsymbol{B}}}_{z}}({k}_{\perp })/{{\rm{P}}}_{\rho }({k}_{\perp })$ ratio of both the FW (in violet) and the A/KAW (orange) power spectra, together with the ratio of the spectra integrated over all the frequencies (in black). The corresponding theoretical ratios (see Appendix E), calculated using the exact wave solutions (dashed curves) of the linearized Hall-MHD equations (Pucci et al. 2016) nicely match both the FW and the A/KAW ratios from the simulation, further confirming the wave nature of the ridges. The approximate formula ${{\rm{P}}}_{{{\boldsymbol{B}}}_{z}}({k}_{\perp })/{{\rm{P}}}_{\rho }({k}_{\perp })={\left({\rm{\Gamma }}\beta \right)}^{2}/4$ for the KAW ratios (see Appendix D) matches only partially both the exact and the computed KAW ratios, due to the fact that such formula is only valid for large angles θ_Bk ≳ 88°. Interestingly, the approximate KAW ratio and the slow/ion-cyclotron (S/IC) exact ratio match well with the ratio of the total spectra at ion-kinetic scales, thus suggesting the presence of almost-perpendicular low-frequency KAWs and/or S/IC. To verify such a possibility, we calculated the ion-velocity polarization ratios (Figure 4(b)). Again, the ratio from the A/KAW ridge nicely match the exact relation at ion-kinetic scales. Moreover, at MHD scales, the ratio of the total spectra (black solid curve) follows the exact A/KAW ratio, as expected since turbulence is Alfvénic at those scales. At ion-kinetic scales, however, neither the approximate KAW (black-dotted line) nor the exact S/IC (blue-dashed curve) ratio match with the computed ratio (solid black curve). This rules out the presence of (energetically relevant) low-frequency KAWs and S/IC and point to a different origin for the matching of the approximate KAW ratio with the ratio of the total spectra in Figure 4(a) (e.g., a state of pressure balance at ion-kinetic scales, see Appendix D).

Zoom In Zoom Out Reset image size

Our results are constrained by the 2D geometry. Wave properties are not completely captured in a 2D model, due to the limited access to the parallel wavenumber (k_∣∣) space. However, we expect that our results are to some extent relevant to a more realistic full 3D case. The spectral turbulent properties measured in 3D simulations are well reproduced by 2D ones (Franci et al. 2018). Moreover, there are regions in our 2D domain where the local mean magnetic field has an in-plane component, thus allowing for the k_∣∣-space to be partially accessible, especially at small scales. This is confirmed by the existence of A/KAW and FW ridges satisfying dispersion relations with θ_Bk < 90° (see Figure 1). Finally and more importantly, theoretical models of KAW turbulence, which require large angles of propagation (k_∣∣ ≪ k_⊥), would be captured by our 2D simulation. In any case, to confirm our findings, we need to extend this study to a fully 3D numerical setup.

To summarize, our results show that 2D Hall-MHD turbulence at sub-ion scales is shaped by low-frequency features (as reported by Parashar et al. 2010). This strongly points to a scenario where the plasma dynamics is dominated by highly intermittent slowly evolving structures where enhanced dissipation occurs (Yang et al. 2017; Camporeale et al. 2018; Papini et al. 2020). Our results indicate that 2D Hall-MHD turbulence on sub-ion scales is not well described in terms of KAW turbulence phenomenology (Howes et al. 2008; Schekochihin et al. 2009; Boldyrev et al. 2013). KAW ratios (Chen et al. 2013a) are not a good proxy for discriminating KAW activity, because they only approximate the properties of KAWs and are also sensitive to the presence of S/IC waves and to pressure-balanced structures (see Appendix D and Verscharen et al. 2017; Grošelj et al. 2019).

Fundamental processes, such as wave–particle resonances and interactions (Markovskii et al. 2006; Sahraoui et al. 2010; Sulem et al. 2016), as well as nongyrotropic effects (Del Sarto et al. 2016), may also contribute to dissipation and heating in space plasmas. Such processes are not properly captured by Hall-MHD models. In this context, future studies extending the FTFIF analysis to large 3D simulations that retain a fully kinetic description of plasmas (Wan et al. 2015; Franci et al. 2018; Grošelj et al. 2019; Franci et al. 2020) will provide fundamental insights into the dynamics and heating in turbulence.

M.P. thanks the Italian Space Agency for the financial support under the contract ASI "LIMADOU scienza +" n° 2020-31-HH.0. A.C. is a member of the Italian "Gruppo Nazionale di Calcolo Scientifico" (GNCS) of the Istituto Nazionale di Alta Matematica "Francesco Severi" (INdAM). L.F. is supported by the UK Science and Technology Facilities Council (STFC) grant ST/T00018X/1. This work was supported by the Programme National PNST of CNRS/INSU cofunded by CNES. We acknowledge partial funding by "Fondazione Cassa di Risparmio di Firenze" under the project HIPERCRHEL. We acknowledge the computing center of CINECA and INAF, under the coordination of the "Accordo Quadro MoU per lo svolgimento di attività congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione", for the availability of computing resources and support (project INA20_C6A55) E. Papini acknowledges CINECA for awarding access to HPC resources under the ISCRA initiative (grant HP10C2EARF and HP10C4A4M2). The MATLAB code implementing the FIF algorithm is available at www.cicone.com. The authors thank L. Matteini and V. Montagud-Camps for useful discussion.

Iterative Filtering (IF; Lin et al. 2009; Cicone et al. 2016) is a technique for the analysis of nonlinear nonstationary signals. IF decomposes a given L² signal f(t) into N simple oscillating functions ${\widehat{f}}_{m}(t)$ called IMCs, such that

$$\begin{eqnarray}&&f(t)=\sum _{m=1}^{N}{\widehat{f}}_{m}(t)+r(t),\end{eqnarray} \tag{ A1 }$$

where r(t) is the residual (or trend) of the decomposition. IF uses a low-pass filter, applied iteratively, to extract the moving average of the signal at a given timescale τ_m , in order to isolate a fluctuating component whose average frequency ω_m ∼ 2π/τ_m is well behaved. Each IMC is given by

$$\begin{eqnarray}&&{\widehat{f}}_{m}(t)=\mathop{\mathrm{lim}}\limits_{n\to \infty }{{ \mathcal S }}_{{\tau }_{m}}^{n}\left[f(t)-\sum _{l=1}^{m-1}{\widehat{f}}_{l}(t)\right],\end{eqnarray} \tag{ A2 }$$

where ${{ \mathcal S }}_{{\tau }_{m}}={ \mathcal I }-{{ \mathcal L }}_{{\tau }_{m}}$, ${ \mathcal I }$ is the identity operator, ${{ \mathcal L }}_{{\tau }_{m}}$ is the integral operator associated to the low-pass filter (Papini et al. 2020), and where n denotes the number of times the operator ${S}_{{\tau }_{m}}$ is applied. τ_m is different for each IMC and increasing with m. Therefore, IMCs with increasing m will contain smaller frequencies ω_m . The low-pass filter operator reads

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\tau }_{m}}[s(t)]={\int }_{-{\tau }_{m}}^{{\tau }_{m}}s(t+\tau ){w}_{m}(\tau )d\tau ,\end{eqnarray} \tag{ A3 }$$

where w_m (τ) [−τ_m , τ_m ] is the kernel function associated to the filter, with compact support [−τ_m , τ_m ]. Here we employ a Fokker–Plank filter (Cicone et al. 2016). This approach has been recently accelerated in what is known as FIF (Cicone 2020).

The procedure used to calculate the k ω-diagram is the following. We consider a field F(x, y, t), with (x,y) defined in a 2D periodic domain on a discrete grid of N_x × N_y points and with t ∈ [t₀, t₁]. F(x, y, t) can be either a scalar field (e.g., density) or a vector component (e.g., B_z ). The spacetime analysis, employs a mixed approach. We perform an FT in space and we then apply FIF in time. The resulting FTFIF decomposition consists of the following steps:

• 1.  
  Perform an FT of F(x, y, t) to obtain its complex Fourier spectrum $\widetilde{F}({k}_{x}^{i},{k}_{y}^{j},t)={\widetilde{F}}^{{ij}}(t)$, with i(j) = 1, ...N_x (N_y ).
• 2.  
  For each real (imaginary) component F^ij (t) of ${\widetilde{F}}^{{ij}}(t)$, perform an FIF decomposition to obtain N^ij IMCs:

  $$\begin{eqnarray}&&{F}^{{ij}}(t)=\sum _{m=1}^{{N}^{{ij}}}{\widehat{F}}_{m}^{{ij}}(t)+{r}^{{ij}}(t).\end{eqnarray} \tag{ B1 }$$
• 3.  
  For each pair (i,j), increase the orthogonality of the set $\{{\widehat{F}}_{m}^{{ij}}(t)\}$ by calculating the cross-correlation matrix

  $$\begin{eqnarray*}&&{\mathrm{Corr}}_{l,n}=\langle {\widehat{F}}_{l}^{{ij}}\cdot {\widehat{F}}_{n}^{{ij}}\rangle /\sqrt{\langle {\left({\widehat{F}}_{l}^{{ij}}\right)}^{2}\rangle \langle {\left({\widehat{F}}_{n}^{{ij}}\right)}^{2}\rangle },\end{eqnarray*}$$

  where

  $$\begin{eqnarray*}&&\langle {\widehat{F}}_{l}^{{ij}}\cdot {\widehat{F}}_{n}^{{ij}}\rangle ={\int }_{{t}_{0}}^{{t}_{1}}{\widehat{F}}_{l}^{{ij}}(t){\widehat{F}}_{n}^{{ij}}(t){dt}.\end{eqnarray*}$$

  Sum the IMCs for which Corr_l,n > 0.6, to obtain a reduced set $\{{\widehat{F}}_{m}^{{ij}}(t)\}$ (with m = 1, M^ij ), which is almost orthogonal.
• 4.  
  For each IMC, calculate its average frequency ${\omega }_{m}^{{ij}}$ and its amplitude ${A}_{m}^{{ij}}=\sqrt{\langle {\left({\widehat{F}}_{m}^{{ij}}\right)}^{2}\rangle }$. At the end you have a set of frequencies and amplitudes $\{{\omega }_{m}^{{ij}},{A}_{m}^{{ij}}\}$, with m = 1, M^ij and i(j) = 1, ...N_x (N_y ).
• 5.  
  Define an (arbitrary) equidistant k ω-grid and interpolate the energy ${\left({A}_{m}^{{ij}}\right)}^{2}$ of each IMC to the four points nearest to $\left({k}_{\perp }^{{ij}}=\sqrt{{k}_{x}^{i,2}+{k}_{y}^{j,2}},{\omega }_{m}^{{ij}}\right)$, such that energy is conserved:

  $$\begin{eqnarray*}&&{\left({A}_{m}^{{ij}}\right)}^{2}=\sum _{{k}_{\perp }={k}_{1},{k}_{2}}\sum _{\omega ={\omega }_{1},{\omega }_{2}}{ \mathcal B }\left[{\left({A}_{m}^{{ij}}\right)}^{2}\right]({k}_{\perp },\omega ),\end{eqnarray*}$$

  where ${ \mathcal B }$ is the interpolation operator and all pair combinations (k_1,2, ω_1,2) denote the coordinates of the four nearest points.
• 6.  
  Finally, sum over i, j, and m to obtain the k ω-power spectrum of F

  $$\begin{eqnarray}&&\begin{array}{l}{{\rm{P}}}_{F}({k}_{\perp },\omega )\\ \quad =\,\displaystyle \frac{1}{{\rm{\Delta }}T}\displaystyle \sum _{i=1}^{{N}_{x}}\displaystyle \sum _{j=1}^{{N}_{y}}\displaystyle \sum _{m=1}^{{M}^{{ij}}}{ \mathcal B }\left[{\left({A}_{m}^{{ij}}\right)}^{2}\right]({k}_{\perp },\omega ),\end{array}\end{eqnarray} \tag{ B2 }$$

  where ΔT = t₁ − t₀ is the length of the temporal interval considered.

The pseudospectral code employed in this work (Papini et al.2019) solves the fully compressible viscous-resistive Hall-MHD equations (in a nondimensionalized form)

$$\begin{eqnarray}&&{\partial }_{t}\rho =-{\rm{\nabla }}\cdot (\rho {\boldsymbol{u}}),\end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\rho {{\rm{d}}}_{t}{\boldsymbol{u}} & = & -{\rm{\nabla }}P+({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}\\ & & +\mu \left[{{\rm{\nabla }}}^{2}{\boldsymbol{u}}+{\rm{\nabla }}({\rm{\nabla }}\cdot {\boldsymbol{u}})/3\right],\end{array}\end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{d}}}_{t}T=\,({\rm{\Gamma }}-1)\left\{-({\rm{\nabla }}\cdot {\boldsymbol{u}})T+\eta | {\rm{\nabla }}\times {\boldsymbol{B}}{| }^{2}/\rho \right.\\ \quad \left.\,+\mu {\rho }^{-1}\left[{\left({\rm{\nabla }}\times {\boldsymbol{u}}\right)}^{2}+4/3{\left({\rm{\nabla }}\cdot {\boldsymbol{u}}\right)}^{2}\right]\right\},\end{array}\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}{\boldsymbol{B}}=\,{\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}\right)+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}}+\\ \,-{d}_{i}/L{\rm{\nabla }}\times [({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}/\rho ],\end{array}\end{eqnarray} \tag{ C4 }$$

where d_t = ∂_t + u · ∇, Γ = 5/3 is the adiabatic index and {ρ, u , B , T, P} denote density, ion-fluid velocity, magnetic field, temperature, and pressure, respectively. Thermodynamic variables are related through the equation of state P = ρ T. All quantities are normalized with respect to a characteristic length L (set to be equal to the ion inertial length d_i ), a plasma density ρ₀, a magnetic field amplitude B₀, the Alfvén velocity ${c}_{A}={B}_{0}/\sqrt{4\pi {\rho }_{0}}={{\rm{\Omega }}}_{i}{d}_{i}$, a pressure ${P}_{0}={\rho }_{0}{c}_{A}^{2}$, and a plasma temperature T₀ = (k_B /m_i )P₀/ρ₀. Ω_i = eB₀/m_i c is the ion-cyclotron angular frequency and m_i is the ion's mass. With this normalization, the Alfvén time is ${\tau }_{A}={d}_{i}/{c}_{A}={{\rm{\Omega }}}_{i}^{-1}$, dynamic viscosity and magnetic resistivity are in units of d_i c_A ρ₀ and d_i c_A , respectively, and the Hall coefficient d_i /L = 1 in Equation (C4). In the above equations, the magnetic field B is expressed in units of the Alfvén velocity. In code units, the electron fluid velocity is u _e = u _i − ∇ × B /ρ.

Equations (C1)–(C4) are solved in a 2.5D periodic Cartesian (x,y) domain, by employing FT for dealiasing and calculating spatial derivatives and by using a third-order Runge–Kutta method for time integration. In 2.5D codes, the space of the coordinates is 2D, vector quantities retain all three xyz-components, and ∂_z = 0 for all variables. The simulation analyzed in the present work employs a 2D box of size L_x × L_y = 128d_i × 128d_i , with 1024² gridpoints and a resolution Δx = Δy = d_i /8. Finally, we set the plasma β = 2P₀/B₀ = 2 and μ = η = 0.001.

Following Boldyrev et al. (2013), the ratio ${\left(\delta {B}_{z}(k)/\delta \rho (k)\right)}^{2}$ for KAWs, obtained in the limit k_z ≪ k_⊥, is given by the formula

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {B}_{z}(k)/{B}_{0}}{\delta \rho (k)/{\rho }_{0}}\right)}^{2}={\left[\displaystyle \frac{{v}_{{T}_{i}}^{2}}{{c}_{A}^{2}}\left(1+\displaystyle \frac{{T}_{e}}{{T}_{i}}\right)\right]}^{2},\end{eqnarray} \tag{ D1 }$$

where T_i (T_e ) is the ion (electron) temperature, and ${v}_{{T}_{i}}$ is the thermal speed of the ions. The above expression can be adapted to the case of adiabatic HMHD by substituting the ion thermal speed with the adiabatic sound speed (${v}_{{T}_{i}}\to {c}_{s}$) and by using the definition $\beta =(2/{\rm{\Gamma }}){c}_{s}^{2}/{c}_{A}^{2}$ together with T_e /T_i = β_e /β_i (where we assumed charge neutrality n_e = n_i = ρ/m_i ) to give

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {B}_{z}(k)/{B}_{0}}{\delta \rho (k)/{\rho }_{0}}\right)}^{2}={\left[\displaystyle \frac{{\rm{\Gamma }}}{2}({\beta }_{i}+{\beta }_{e})\right]}^{2}=\displaystyle \frac{{\left({\rm{\Gamma }}\beta \right)}^{2}}{4},\end{eqnarray} \tag{ D2 }$$

where $\beta ={\beta }_{i}+{\beta }_{e}=(2/{\rm{\Gamma }}){c}_{s}^{2}/{c}_{A}^{2}$ and ${c}_{s}^{2}={\rm{\Gamma }}{P}_{0}/{\rho }_{0}$. We note that the above formula can be obtained self-consistently from the KAW linear solutions of ideal HMHD (see Equation (E2)) in the limit k_z ≪ k_⊥.

Equation (D2) is also fulfilled by slowly evolving pressure-balanced structures. In Fourier k ω-space, in the limit of low frequencies and for wavenumbers smaller than the inverse Kolmogorov dissipation scale, Equation (C2) reads

$$\begin{eqnarray}&&\delta P(k)+\displaystyle \frac{\delta {B}_{z}(k){B}_{0}}{4\pi }\simeq 0\end{eqnarray} \tag{ D3 }$$

to first order in the turbulent fluctuations. By using $\delta P={c}_{s}^{2}\delta \rho$, we then obtain

$$\begin{eqnarray}&&{c}_{s}^{2}\delta \rho (k)\simeq -\displaystyle \frac{\delta {B}_{z}(k){B}_{0}}{4\pi }.\end{eqnarray} \tag{ D4 }$$

The latter equation may be rearranged to give the ratio between the parallel magnetic fluctuations δ B_z (k) and the density fluctuations δ ρ(k), to obtain

$$\begin{eqnarray}&&\displaystyle \frac{\delta {B}_{z}{\left(k\right)}^{2}}{\delta \rho {\left(k\right)}^{2}}={\left(\displaystyle \frac{{c}_{s}^{2}}{{c}_{A}^{2}}\displaystyle \frac{{B}_{0}}{{\rho }_{0}}\right)}^{2}={\left(\displaystyle \frac{{\rm{\Gamma }}\beta }{2}\displaystyle \frac{{B}_{0}}{{\rho }_{0}}\right)}^{2},\end{eqnarray} \tag{ D5 }$$

which, in adimensional units and by using ${{\rm{P}}}_{X}(k)={\left(\delta X(k)\right)}^{2}/k$, gives the expression

$$\begin{eqnarray}&&\frac{{{\rm{P}}}_{{B}_{z}}}{{{\rm{P}}}_{\rho }}=\frac{{\left({\rm{\Gamma }}\beta \right)}^{2}}{4}\end{eqnarray} \tag{ D6 }$$

that is equal to Equation (D2) for the approximated KAW ratios.

In presence of the Hall term and below ion scales (kd_i ≳ 1), the three MHD wave solutions (fast, slow, and Alfvén) are modified. The fast mode turns into a whistler mode (FW), the Alfvén mode transitions to a kinetic Alfvén mode (A/KAW), and the slow mode becomes a (compressible) ion-cyclotron mode (S/IC) (e.g., Hollweg 2000; Galtier 2016). The polarization ratios for FW, A/KAW, and S/IC waves are calculated from the eigenmode solutions {δ u (ω, k ), δ B (ω, k ), δ ρ(ω, k )} of the linearized Hall-MHD equations. The expression for the dispersion relation, δ u , and δ B can be found in Pucci et al. (2016). Following their convention, here we write the expression for the amplitude of the density

$$\begin{eqnarray}&&\delta \rho (\omega ,{\boldsymbol{k}})=\left|{\rho }_{0}a\displaystyle \frac{({k}_{| | }^{2}{c}_{A}^{2}-{\omega }^{2}){k}_{\perp }{{\rm{\Omega }}}_{i}}{({\omega }^{2}-{k}^{2}{c}_{s}^{2}){k}_{| | }\omega }\right|,\end{eqnarray} \tag{ E1 }$$

where ω is the eigenfrequency, a is the perturbation amplitude, and k_∣∣(⊥) is the parallel (perpendicular) component of the wavevector k with respect to the background magnetic field, thus $\tan ({\theta }_{{Bk}})={k}_{\perp }/{k}_{| | }$. Note that the definition of parallel direction used here differs from that of the numerical simulation, in which the parallel direction is that of the global mean magnetic field (i.e., the z-direction).

The exact formulas for the polarization ratios shown in this work read

$$\begin{eqnarray}&&{\left(\frac{\delta {B}_{| | }{\rho }_{0}}{\delta \rho {B}_{0}}\right)}^{2}=\frac{{\left({\omega }^{2}-{k}^{2}{c}_{s}^{2}\right)}^{2}}{{k}^{4}{c}_{A}^{4}},\end{eqnarray} \tag{ E2 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{\delta {B}_{\perp }{\rho }_{0}}{\delta \rho {B}_{0}}\right)}^{2}=\displaystyle \frac{{\left({\omega }^{2}-{k}^{2}{c}_{s}^{2}\right)}^{2}}{{\left({k}_{\perp }/{k}_{| | }\right)}^{2}}\\ \quad \left[\displaystyle \frac{{\left(\omega /{{\rm{\Omega }}}_{i}\right)}^{2}}{{\left({\omega }^{2}-{k}_{| | }^{2}{c}_{A}^{2}\right)}^{2}}+\displaystyle \frac{1}{{k}^{4}{c}_{A}^{4}}\right],\end{array}\end{eqnarray} \tag{ E3 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{\delta {u}_{\perp }}{\delta {u}_{| | }}\right)}^{2}={\left[\displaystyle \frac{\omega {k}_{| | }{c}_{A}^{2}({\omega }^{2}-{k}^{2}{c}_{s}^{2})}{{{\rm{\Omega }}}_{i}{k}_{\perp }{c}_{s}^{2}({\omega }^{2}-{k}_{| | }^{2}{c}_{A}^{2})}\right]}^{2}\\ \,+{\left[\displaystyle \frac{({\omega }^{2}-{k}_{| | }^{2}{c}_{s}^{2})}{{k}_{| | }{k}_{\perp }{c}_{s}^{2}}\right]}^{2},\end{array}\end{eqnarray} \tag{ E4 }$$

where the ratio for the desired branch (FW, A/KAW, or S/IC) is obtained by plugging in the corresponding eigenfrequency ω (whose analytical expression can be found in Vocks et al. 1999). The analogous ratios from the numerical simulation have been obtained by using $\delta {B}_{\perp }^{2}/k={{\rm{P}}}_{{B}_{\perp }}$, $\delta {{B}_{| | }}^{2}/k={{\rm{P}}}_{{B}_{z}}$, $\delta {u}_{\perp }^{2}/k={{\rm{P}}}_{{u}_{\perp }}$, and $\delta {{u}_{| | }}^{2}/k={{\rm{P}}}_{{u}_{z}}$. Equations (E2) and (E4) give the exact ratios used in Figure 4.

In Figure 5 we report the polarization ratios for ${{\rm{P}}}_{{{\boldsymbol{B}}}_{\perp }}({k}_{\perp })/{{\rm{P}}}_{\rho }({k}_{\perp })$. Again, as shown in Figure 4, the A/KAW ratio nicely matches the exact relations at ion-kinetic scales. Moreover, at MHD scales, the ratio of the total spectra (black solid curve) follows the exact A/KAW ratio, as expected since turbulence is Alfvénic at those scales. At ion-kinetic scales however, this is not the case. The approximate KAW ratio (horizontal dotted lines) matches the ratio of the total spectra for ${{\rm{P}}}_{{{\boldsymbol{B}}}_{\perp }}/{{\rm{P}}}_{\rho }$, (although such a match is not as good as in the case of ${{\rm{P}}}_{{{\boldsymbol{B}}}_{z}}/{{\rm{P}}}_{\rho }$). However, the presence of energetically relevant low-frequency KAWs has already been ruled out by the $\delta {u}_{\perp }^{2}/\delta {u}_{z}^{2}$ ratio. This casts some doubts on the use of ${{\rm{P}}}_{{{\boldsymbol{B}}}_{\perp }}/{{\rm{P}}}_{\rho }$ as a proxy for KAW activity.

Zoom In Zoom Out Reset image size

The approximate relations for the ${{\rm{P}}}_{{{\boldsymbol{B}}}_{\perp }}/{{\rm{P}}}_{\rho }$ and the ${{\rm{P}}}_{{{\boldsymbol{u}}}_{\perp }}/{{\rm{P}}}_{{{\boldsymbol{u}}}_{z}}$ KAW ratios are

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {B}_{\perp }{\rho }_{0}}{\delta \rho {B}_{0}}\right)}^{2}\simeq \displaystyle \frac{{\rm{\Gamma }}\beta }{2}\left(1+\displaystyle \frac{{\rm{\Gamma }}\beta }{2}\right),\end{eqnarray} \tag{ E5 }$$

$$\begin{eqnarray}&&{\left(\frac{\delta {u}_{\perp }}{\delta {u}_{| | }}\right)}^{2}\simeq \frac{2}{{\rm{\Gamma }}\beta }+1.\end{eqnarray} \tag{ E6 }$$