Power Anisotropy, Dispersion Signature and Turbulence Diffusion Region in the 3D Wavenumber Domain of Space Plasma Turbulence

We explore the multi-faceted important features of turbulence (e.g., anisotropy, dispersion, diffusion) in the three-dimensional (3D) wavenumber domain ($k_\parallel$, $k_{\perp,1}$, $k_{\perp,2}$), by employing the k-filtering technique to the high-quality measurements of fields and particles from the MMS multi-spacecraft constellation. We compute the 3D power spectral densities (PSDs) of magnetic and electric fluctuations (marked as $\rm{PSD}(\delta \mathbf{B}(\mathbf{k}))$ and $\rm{PSD}(\delta \mathbf{E}'_{\langle\mathbf{v}_\mathrm{i}\rangle}(\mathbf{k}))$), both of which show a prominent spectral anisotropy in the sub-ion range. We give the first 3D image of the bifurcation between power spectra of the electric and magnetic fluctuations, by calculating the ratio between $\rm{PSD}(\delta \mathbf{E}'_{ \langle\mathbf{v}_\mathrm{i}\rangle}(\mathbf{k}))$ and $\rm{PSD}(\delta \mathbf{B}(\mathbf{k}))$, the distribution of which is related to the non-linear dispersion relation. We also compute the ratio between electric spectra in different reference frames defined by the ion bulk velocity, that is $\mathrm{PSD}(\delta{\mathbf{E}'_{\mathrm{local}\ \mathbf{v}_\mathrm{i}}})/\mathrm{PSD}(\delta{\mathbf{E}'_{ \langle\mathbf{v}_\mathrm{i}\rangle}})$, to visualize the turbulence ion diffusion region (T-IDR) in wavenumber space. The T-IDR has an anisotropy and a preferential direction of wavevectors, which is generally consistent with the plasma wave theory prediction based on the dominance of kinetic Alfv\'en waves (KAW). This work manifests the worth of the k-filtering technique in diagnosing turbulence comprehensively, especially when the electric field is involved.


INTRODUCTION
Plasma turbulence is ubiquitous in most space environments, including the solar wind (SW) (Alexandrova et al. 2013;Bruno & Carbone 2013) and the magnetosheath (Sahraoui et al. 2003;Huang et al. 2014;Sahraoui et al. 2020). The investigation of turbulence is necessary to understand the acceleration, heating and transports of space plasmas (Schekochihin et al. 2009;Howes et al. 2011).
Turbulence exhibits an anisotropic distribution of fluctuation energy in wavenumber (k) space, i.e., a (power) spectral anisotropy (Horbury et al. 2012;Oughton et al. 2015;Narita 2018). This kind of anisotropy is typically characterized with weaker and stronger power levels when sampled along the parallel and perpendicular directions relative to the background magnetic field, respectively. In earlier studies, the "Maltese-cross" pattern of the magnetic correlation function at 1 au shows that turbulence exhibits features of a "Slab+2D" configuration (Matthaeus et al. 1990), which consists of parallel waves and perpendicular structures. As an alternative, the Critical Balance (CB) theory proposed by Goldreich & Sridhar (1995) assumes that the magnitudes of the linear (e.g. Alfvén) time scale τ A and of the non-linear time scale τ NL are comparable, leading to a balance between the effects of wave propagation and non-linear interaction. Applying the CB theory to Alfvén waves yields scaling indices of power spectral density (PSD): -2 in the parallel direction and -5/3 in the perpendicular direction, which is in agreement with the observed dependence of PSD indices on the sampling direc-tion relative to the local background magnetic field (θ V B ) in single-point timeseries analyses of the the SW (Horbury et al. 2008). The observation of the scaling anisotropy at Magnetohydrodynamics (MHD) scales is to some degree sensitive to the determination of the background magnetic field direction. However, the power anisotropy still exists throughout the inertial range .
Anisotropy also exists in the kinetic range (Leamon et al. 1998;Sahraoui et al. 2010b). Involving appropriate kinetic wave modes in critically balanced cascade models permits the prediction of the anisotropy at kinetic scales. Ideally, a critically balanced Kinetic Alfvén Wave (KAW) cascade would Cho & Lazarian 2004;Chen et al. 2010). Spectral anisotropy of turbulence is also found in the magnetosheath (Alexandrova et al. 2008;He et al. 2011c). A recent statistical and robust survey conducted by Wang et al. (2020) has found a scale-dependent 3D anisotropy pattern down to the electron scales, with the relation l ∝ l 0.72 ⊥ for the parallel and the perpendicular correlation lengths with respect to the background magnetic field in the sub-ion range without being significantly affected by intermittency as previously thought. The observational scaling cannot be explained by the theory based on CB-type KAW, which predicts l ∝ l 1/3 ⊥ (equivalent to k ∝ k 1/3 ⊥ ), and calls for further investigation. The coexistence of different wave modes at kinetic scales (references) complicates the study of small-scale plasma turbulence furthermore (He et al. 2011a,b;Podesta 2012;He et al. 2015;Zhu et al. 2019).
Nevertheless, when dealing with data offered by only one spacecraft, as the spatial and temporal variations are entangled, it is impossible to define the full four-dimensional PSD directly. This problem is a key driver for multi-spacecraft missions like Cluster II (Escoubet et al. 2001) and Magnetospheric Multiscale Mission (MMS) (Burch et al. 2016), and advanced multi-spacecraft analysis techniques, including the k-filtering technique (also called the wave telescope technique) we apply in this work. It was proposed by Capon (1969), and introduced to space physics by Pinçon & Lefeuvre (1991). As a generalized minimum variance analysis, it gives an estimation of the four-dimensional PSD(ω, k) based on Fourier Transformation so that time periodicity and spatial periodicity are separated. Cluster II is the first mission that has offered applicable datasets for the technique (Glassmeier et al. 2001;Sahraoui et al. 2003). Sahraoui et al. (2006) investigated an event dominated by mirror-mode turbulence with strong anisotropy, starting the investigation into turbulence with this technique. Narita et al. (2010b) presented the first 4D PSD of turbulence, which paved the way for our methodology presented in this article. In this method, the smallest separation between the spacecraft defines the smallest scales of the turbulence that can be resolved. The MMS mission provides high-resolution datasets and smaller satellite separation down to the electron scale for this technique. Narita et al. (2016) presented a comprehensive case study of turbulence based on MMS observations, including wave mode, dispersion relation and propagating direction of waves. Roberts et al. (2019) extended the scope of the measurement beyond the magnetic field, e.g. velocity, density and thermal velocity of ions, and estimated power index anisotropy. However, k-filtering analyses of the turbulence electric field are rare. As an early and inspiring work, Tjulin et al. (2005) argued that in order to use more information and to estimate polarization, we must include the electric field in this technique. After Tjulin, however, we can hardly find applications involving the electric field.
The fluctuating electric field and its PSD play important roles in plasma turbulence. According to Faraday's law, the bifurcation of the PSDs of electric fluctuations and magnetic fluctuations is a signature of dispersion at kinetic scales, which was used by Bale et al. (2005) to support that turbulence has a KAW nature at the sub-ion scale. The electric field is often evoked with the current to estimate the energy-transfer rate between particles and fields as J · E or J · E (Zenitani et al. 2011;He et al. 2019;Duan et al. 2020). Additionally, by computing the electric field in the ion-flow reference frame E = E + v i × B and using |δE |/|δE| as a metric, where δE denotes the fluctuation of the electric field, we determine the plasma demagnetization in magnetic reconnection events (Hesse et al. 1999;Birn & Priest 2007;Lu et al. 2010). This measurement has been adopted by Duan et al. (2018) for the computation of energy transfer in plasma waves based on linear kinetic theory to illustrate the diffusion of magnetic flux relative to the flow of plasma in wavenumber space, which becomes significant when approaching kinetic scales. He et al. (2020) further proposed the concept of turbulence ion and electron diffusion regions (T-IDR and T-EDR) based on the observed the ratio PSD(δE i/e,local )/PSD(δE i,global ) as a reasonable alternative to |δE |/|δE|. In conclusion, the spectrum of the fluctuating electric field is related to basic features and mechanisms of plasma turbulence, especially at sub-ion scales. Reliable electric PSDs are necessary for the investigations of this topic.
Out of the considerations above, we apply the k-filtering technique to a case study of plasma turbulence in the magnetosheath, using measurements from MMS. We compute PSDs from the highresolution magnetic and electric field signals, i.e., PSD(δB(k)) and PSD(δE v i (k)), respectively. We compare these two PSDs to attain the dispersion signature in 3D wavenumber space. We compute the PSD of the electric field in the local ion-flow frame, PSD(δE local v i ), and in the global (averaged) ion-flow frame, PSD(δE v i ), respectively. By calculating the ratio PSD(δE local v i )/PSD(δE v i ), we identify the T-IDR in 3D wavenumber space. In Section 2, We introduce our dataset and methodology. We present our PSDs of the magnetic and electric fluctuations, and discuss the dispersion signature in Section 3. The T-IDR will be shown and discussed in Section 4.

The K-filtering Technique and Dealing with the Doppler Effect
Here we briefly review the k-filtering technique (Paschmann & Schwartz 2000). We use the original form of the technique, which consists of three main steps. Firstly we apply a Fourier Transformation to the signals, and then combine these resulting Fourier images into a vector A(ω, r i ), e.g. when we choose the magnetic field signals: where the signal A can be a combination of vector components (E, B or [B, E] T ), or scalars (|B|, n, T , etc.). Then we combine the A(ω, r i ) from the four different satellites together to [A(ω, r 1 ), A(ω, r 2 ), A(ω, r 3 ), A(ω, r 4 )] T (abbreviated as A(ω) with r i omitted), and with proper averaging, we have the covariance matrix M A (ω), We introduce the steering (wave-propagating) matrix H(k), where I denotes the unit matrix, Using the Lagrange multiplier technique, we obtain for the estimation of PSD the following expression, Previous studies remind us of several limitations of this technique ).The technique is based on the assumption of quasi-stationarity and quasi-homogeneity of the turbulent fluctuations. Moreover, the geometric configuration of the spacecraft constellation determines the first Brillouin zone where the estimation is reliable. Due to the spatial aliasing effect that fluctuations separated by ∆k = 2nπ/d in k-space cannot be distinguished. The maximum wavenumber k max resolved in this technique is π/d, where the d is the appropriate (averaged) separation of spacecraft in the constellation. This wavenumber k max corresponds to a wavelength of λ = 2d. A configuration of the constellation close to a regular tetrahedron maximizes the Brillouin zone volume and reduces angular aliasing, which would otherwise create fake anisotropy (Narita 2009). Empirically given by Sahraoui et al. (2010a) to guarantee accurate positions of power peaks, the minimum resolved wavenumber k min corresponds to ∼ π/5d.
Due to the Doppler effect, there is a shift between the frequency of fluctuations measured in the spacecraft frame (ω sc ) and the frequency in the reference frame defined by the plasma flow (ω pl ).
In order to compare with theoretical predictions and eliminate the Doppler shift, we reconstruct PSD(ω pl , k) from the direct result PSD(ω sc , k). Using the Doppler relation ω pl = ω sc − k · v, we map the value of PSD(ω sc , k) directly to the corresponding PSD(ω pl , k). As the relation ω pl = ω sc − k · v corresponds to oblique planes in the 4D (ω sc , k) space if ω pl is set, the reconstruction method can be symbolized like cutting 4D bread into pieces with an oblique knife. In the computation, independent variables are discrete, so linear interpolation is employed here. When ω pl is linked with negative ω sc , the equality PSD(ω, k) = PSD(−ω, −k) is used. A figuratively description of this procedure is included in the Appendix.

Dataset for Analysis
The observation analysed is from MMS from 9:24:11 to 9:25:07, October 16th, 2015, when the satellites are located in the magnetosheath close to the dusk magnetopause, at a distance of 11.9R E from the earth. The same event has been reported by Chen & Boldyrev (2017), as an event with anisotropy and kinetic Alfvén nature in the sub-ion range, but we neglect the last 17 seconds because the background magnetic field B 0 rotated by about 10 degrees then. The magnetic field data is from the FGM instrument (Russell et al. 2016). The electric field is measured by the SDP (Lindqvist et al. 2016) and ADP (Ergun et al. 2016) instruments. We use data from the FPI instrument for density, bulk velocity and temperature of the protons (Pollock et al. 2016). The timeseries of the event is shown in Figure 1, with all the variables in the Geocentric Solar Ecliptic (GSE) coordinate system. Some key parameters are listed in Table 1. The Alfvén speed exceeds the average ion flow velocity, which potentially causes a violation of the assumptions underlying Taylor's hypothesis. The k-filtering technique does not depend on Taylor's hypothesis, though.

SIGNATURE
The analysed time interval fulfills the conditions for the application of the k-filtering method. The satellites almost formed a regular tetrahedron with Elongation = 0.09 and P lanarity = 0.17 (Robert et al. 1998). The average separation of the satellites, d avg , sets a k max ∼ 0.21km −1 , consistent with k max d i ∼ 12, and k min ∼ 0.04km −1 , consistent with k min d i ∼ 2.4. We rotate the reference frame so that B 0 = B 0ê , and setê ⊥1 andê ⊥2 along x GSE × B 0 and B 0 ×ê ⊥1 , respectively. Figure 2(a-f) show the reduced 2D PSDs of the magnetic and the electric fields, PSD(δB) and PSD(δE v i ), abbreviated as P B and P E v i henceforth. To reduce the 4D spectrum to 2D in the (k i , k j ) plane, we integrate the 4D PSD over ω pl and the other component of k. In these panels, we find a spectral anisotropy in the wavenumber space for both the magnetic and electric fluctuations, although the electric fluctuations are less anisotropic. Figure 2 gives the first 3D image of the bifurcation between P B and P E v i .

Panels (a) and (d)
give an indication of a non-axisymmetry in the power distribution, showing elongated power distributions. This pattern of P (k ⊥1 , k ⊥2 ) reflects that in this event with a relatively short time interval, the kinetic waves have limited directions of k ⊥ , which do not distribute evenly in the whole range of azimuthal angle. We compute the P E v i /P B ratio and normalize it to the Alfvén speed, which is presented in Figure 2(g-i). The ratio is well below the unity at small k (large scales), and increases to values greater than unity when |k| increases along the direction perpendicular to B 0 , where the fluctuation energy concentrates. The ratio increases even more steeply along the directions not perpendicular to B 0 . It exceeds 10 at |k| = 10 when k k ⊥ .
In Figure 3, we present P B , P E v i and 1 V 2 in 3D wavenumber space. The reduced PSDs in wavenumber space are shown after integration over ω pl . All properties mentioned in the above 2D plots are also visible in these 3D plots. We further compute the ratio of P E local v i to P E v i to measure the demagnetization and the diffusion of ions. These two PSDs are computed from electric fields in the local ion-flow frame, where the angle brackets mean averaging over the whole interval, respectively. Figure 4 shows the ratio after appropriate integrations.
As the figure shows, the ratio is lower in regions close to quasi-perpendicular wavevectors (k ⊥ ) than near quasi-parallel wavevectors (k ). In panel (a), the region with a low ratio (less than 1) covers a large range of the parameter space. In panel (b), the low ratio range reaches k ⊥2 d i = 5, and unlike the PSD, it is far less symmetric in k ⊥2 . The corresponding 3D version is shown in Figure 5.  Figure 4.
where P denotes PSD integrated over ω pl and different components of k as defined in Figure 2.
(a) (b) Figure 5. Reduced 3D distribution of the ratio P E local v i /P E v i after integration over ω pl .
Using the New Hampshire Dispersion Relation Solver (NHDS) (Verscharen & Chandran 2018) based on linear Vlasov-Maxwell theory, we calculate the fluctuations δB, δE, δv i and δv e of eigenmodes belonging to the Alfvén-mode branch under the same plasma condition as measured (see Table   1). Figure 6 is derived from the output of NHDS. Figure 6(a) shows the modelled distribution of |δE local v i |/|δE| in wavenumber space in the plasma frame, which is considered equivalent to the square The ratio |δE local v i |/|δE| stays below 1.0 when k d i < 0.5 and k ⊥ d i < 1.0. It exceeds unity at smaller scales. A sharp edge near θ kB = 30 • is visible in panel (a), which separates the high and low ratios at small and large θ kB , respectively. Another ridge-like boundary between low and high ratios is also found at θ kB = 70 • when kd i 2.5. Figure 6(b) displays the normalized magnetic helicity σ m . We also see an edge at θ kB ∼ 30 • in this panel corresponding to σ m = 0. The normalized magnetic helicity σ m is negative (positive) at θ kB < 30 • (θ kB > 30 • ). The positive σ m at large θ kB is an evident feature of KAW, which is distinguished from the negative σ m at small θ kB characteristic for ion cyclotron waves (ICWs).
According to Figure 6(a), when ICWs dominate, the metric of diffusion, 1.0 significantly. For the ICWs at ion scales, the fluctuating electric field δE can be approximated as −δv e × B 0 , since the electrons are still frozen-in with the magnetic fields. Therefore, the electric field in the local ion-flow frame can be approximated as δE ∼ (−δv e + δv i ) × B 0 . On the other hand, we know that |δv i |>|δv e | and |δv i − δv e |>|δv e | for ICWs, which means that ions are the primary current carrier of the wave-current density, according to the plasma wave theory. Therefore, we anticipate that |δE |>|δE| for the ICWs. Comparing to this theoretical prediction, we observe that the k-filtering results do not show such a high ratio of the ion diffusion metric at quasi-parallel wavenumbers. This finding suggests that ICWs are not detectable in our measurement interval.

SUMMARY AND DISCUSSION
For a case study of turbulence in the terrestrial magnetosheath, we compute the 4D PSDs of the turbulent magnetic and electric fields, recalculate the frequency in the plasma frame from its counterpart in the spacecraft frame, and reconstruct the PSDs with the frequency in the plasma frame. Both P B and P E v i show the typical anisotropy in the 2D plane of (k , k ⊥ ), but P E v i is less anisotropic. We also find some degree of non-axisymmetry in the (k ⊥1 , k ⊥2 ) plane, which agrees with Roberts et al. (2019). Such non-axisymmetry with some elongation of the PSD suggests that oblique kinetic Alfvén waves at various k ⊥ are concentrated within a finite range of angles instead of the 2π omni-direction.
The bifurcation between P B and P E v i was evaluated through For the first time, we computed the ratio between PSD(δE local v i ) in the local ion bulk flow frame and PSD(δE v i ) in the global ion bulk flow frame. We computed the electric field power ratio to identify the ion demagnetization effect in wavenumber space at kinetic scales and compared it with the prediction of ion demagnetization in the Alfvén wave branch as calculated by NHDS. The ratio of the event exhibits anisotropy and asymmetry of wavevectors. Considering the strengths of this method as described above, we anticipate that the k-filtering technique will become a powerful technique for the comprehensive investigation of diffusion physics in turbulence.
Although the k-filtering technique has shown considerable capability, it suffers from limitations and errors, as mentioned above. Especially spatial aliasing is a central issue, the challenge of which has been theoretically discussed by Narita (2009). The spatial aliasing effect prevents us from quantifying the distribution of PSDs in (ω, k) and quantitatively comparing reduced PSDs with those obtained directly from Fourier/wavelet transformations of the time sequences. The limitation of this method to around (or even less than) one decade in wavevector space is another weakness of this technique, which prevents the simultaneous analysis of different plasma scales. A more flexible constellation with more than 4 satellites would be a prospect to be anticipated (Zhang et al. 2020;Dai et al. 2020).
Unfortunately, phase spectra of variables in (ω, k) space cannot be acquired through the k-filtering technique. As a consequence, we can not further study the power spectral distribution of different component disturbances (e.g., δE and δE ⊥ ), nor can we study the combined spectra of different disturbances (e.g., δJ · δE for the dissipation rate spectrum).
We sincerely thank Olga Alexandrova for the discussion with us at the vEGU meeting 2021.

APPENDIX
We described our method with Figure 7, which imitating the presentation of Figure 9 in the paper Narita et al. (2010a). Panel(a) is a slice cut from the PSD B (ω sc , k) at k x = 0.0043km −1 and k z = 0.0087km −1 . After the Doppler shift correction, the slice's shape changes from a rectangle to a parallelogram, which is shown in panel(b). For a proper integration following, we require a box-shaped PSD B (ω pl , k), so we reserve the distribution in the red rectangle in panel(b) and append the grey patch from its mirror counterpart using the equality mentioned above. Panel(c) shows the resulted slice with frequencies in the plasma frame. The 4D PSD B (ω pl , k) is the ensemble of all the slices like this one.