Observational quantification of three-dimensional anisotropies and scalings of space plasma turbulence at kinetic scales

A statistical survey of spectral anisotropy of space plasma turbulence is performed using five years measurements from MMS in the magnetosheath. By measuring the five-point second-order structure functions of the magnetic field, we have for the first time quantified the three-dimensional anisotropies and scalings at sub-ion-scales ($<$ 100 km). In the local reference frame $(\hat L_{\perp}, \hat l_{\perp}, \hat l_{\parallel})$ defined with respect to local mean magnetic field $\boldsymbol {B}_0$ (Chen et al. 2012), the"statistical eddies"are found to be mostly elongated along $\boldsymbol {B}_0$ and shortened in the direction perpendicular to both $\boldsymbol {B}_0$ and local field fluctuations. From several $d_i$ (ion inertial length) toward $\sim$ 0.05 $d_i$, the ratio between eddies' parallel and perpendicular lengths features a trend of rise then fall, whereas the anisotropy in the perpendicular plane appears scale-invariant. Specifically, the anisotropy relations for the total magnetic field at 0.1-1.0 $d_i$ are obtained as $l_{\parallel} \simeq 2.44 \cdot l_{\perp}^{0.71}$, and $L_{\perp} \simeq 1.58 \cdot l_{\perp}^{1.08}$, respectively. Our results provide new observational evidence to compare with phenomenological models and numerical simulations, which may help to better understand the nature of kinetic scale turbulence.


INTRODUCTION
The energy distribution at a certain scale (or k space) is known to be not isotropic in the turbulence of magnetized plasma, also known as spectral anisotropy (Cho & Vishniac 2000). This particular feature reflect the preferential direction of the energy cascade with respect to the local background magnetic field B 0 (Podesta 2009). Most of our experimental knowledge of space plasma anisotropy comes from in-situ observations made within the solar wind (SW), which is a nearly collisionless plasmas stream released from the Sun (Bruno & Carbone 2013).
At large magnetohydrodynamic (MHD) scales, the pattern of correlation function for the magnetic field at 1 AU has two major components referred to as "Maltese cross", exhibiting elongations in both parallel and perpendicular direction with regard to B 0 (Matthaeus et al. 1990). This signature is summarized as the "slab+2D" model, which assumes no specific nature of the fluctuations but just describe the fluctuations as a combination of waves with k and structures with k ⊥ . Another type of anisotropy model is based on "critical-balance (CB)" conjecture, where the key hypothesis relies on the comparable scale of the linear Alfvén time and turbulence non-linear time in a vanishing crosshelicity system (Goldreich & Sridhar 1995). As a result, the spectral anisotropy scales as k ∝ k 2/3 ⊥ . By introducing the idea of "dynamic alignment" between magnetic and velocity field fluctuations, Boldyrev (2006)  Indeed, numerous observations have found agreement between measurements and CB theories (Horbury et al. 2008;Luo & Wu 2010;Chen et al. 2011Chen et al. , 2012. Despite this consistency, recent revisit of anisotropy in the solar wind have reported some puzzling results and raised more concerns to be considered, such as intermittency (Wang et al. 2014;Pei et al. 2016;Yang et al. 2017), the discrepancy between velocity field and magnetic field anisotropy (Wicks et al. 2011;Wu et al. 2019a,b;Yan et al. 2016), dependence on the heliocentric distance (He et al. 2013) and solar wind expansion (Verdini et al. , 2019. Moreover, the 3D self-correlation functions are shown to be isotropic in (Wang et al. 2019b;Wu et al. 2019a,b). Therefore, the anisotropic nature of the solar wind at MHD scales is still an open question.
At kinetic scales, the turbulence still remains or become much anisotropic (i.e., Chen et al. (2010); Oughton et al. (2015) and references therein). The standard Kinetic Alfvén Wave (KAW) turbulence model, also on basis of CB conjecture by assuming the linear KAW propagation time to be comparable to nonlinear time, predicts an anisotropy scaling of k ∝ k 1/3 ⊥ (Howes et al. 2008;Schekochihin et al. 2009). It is also important to pay attention to the much complicated physical process at kinetic scales than at MHD regime due to plasma kinetic effects (see the review by Alexandrova et al. (2013) and Alexandrova et al. (2020)). For example, the modified KAW turbulence model with intermittent 2D structures has k ∝ k 2/3 ⊥ (Boldyrev & Perez 2012). Zhao et al. (2016) suggested a model for kinetic-scale Alfvénic turbulence which incorporate the dispersion and intermittency effects. Boldyrev & Loureiro (2019) considered the decisive role played by the tearing instability in setting the aspect ratio of eddies and hence predicted the spectral anisotropy scalings between k k 2/3 ⊥ and k k ⊥ . Most recently, Landi et al. (2019) proposed a phenomenological model considering the intermittent two-dimensional structures in the plane perpendicular to B 0 . In their model, the prescribed perpendicular aspect ratio of these structures could determine the anisotropy as k ∝ k 1/3(α+1) ⊥ , where α is proportional to the space-filling of the turbulence.
In recent high resolution three-dimensional kinetic simulations, the spectral anisotropy has received considerable attentions (Grošelj et al. 2018;Franci et al. 2018;Arzamasskiy et al. 2019;Cerri et al. 2019;Landi et al. 2019). Based on different methods of measuring the anisotropy, dissimilar scaling relations have been found (i.e., k ∝ k Grošelj et al. (2018) and k ∝ k ⊥ in Arzamasskiy et al. (2019). Specifically, for the analysis based on multi-point local structure functions (SF ) which will be introduced in Section 2, the anisotropy tends to become "frozen" when approaching ion scales (i.e., k ∝ k 0.8 ⊥ in Landi et al. (2019)). By comparing the results from three different simulations including the hybrid particle-in-cell (PIC), Eulerian hybrid-Vlasov, and fully kinetic PIC codes, Cerri et al. (2019) found that the anisotropy scalings tend to converge to l ∝ l 2/3 ⊥ based on a unified analysis of five-point SF s. The Earth's magnetosheath (MSH) offer a unique lab different from SW, such as enhanced compressibility, intermittency, as well as the kinetic instabilities (Alexandrova 2008). In addition, the spacecraft measurements in MSH tend to cover a wider range of angle between bulk flow velocity and B 0 as compared with solar wind (He et al. 2011), thus allowing us to diagnose 3D nature of the fluctuations in a relatively short interval. Using Cluster measurements, Mangeney et al. (2006) found the strong anisotropy of the electromagnetic fluctuations, with k ⊥ extending for two decades within the kinetic range kd e ∈ [0.3, 30]. Sahraoui et al. (2006) showed the anisotropic behaviour up to kρ i = 3.5 in a mirror structure event. Alexandrova et al. (2008) surveyed 6 events and found the dominance of 2D turbulence (k ⊥ k ) above the spectral break in the vicinity of ion scale. In addition, due to the Doppler shift, the magnetic fluctuations have more energy along the B × V direction in their analysis. Similar 2D turbulence at kinetic scales was observed in a recent statistical study of magnetic field turbulence in the solar wind . He et al. (2011) computed the spatial correlation functions of both magnetic field and density fluctuations in the 2D (l , l ⊥ ) plane, it is shown that the turbulence close to ion scales is comprised of two populations, where the major component is mostly transverse and the minor one is oblique. Using measurements from Magnetospheric Multiscale mission (MMS) (Burch et al. 2016), Chen & Boldyrev (2017) studied the two-point SF of the magnetic field in the same plane and provided evidence of strong anisotropy at smaller scales (11 < kρ i < 57). Another event recorded by MMS show that parameters such as magnetic field, density, ion velocity, and ion thermal speed all exhibit anisotropy in the spectral index up to kρ i ∼ 1 (Roberts et al. 2019).
Despite these case studies, a comprehensive in-situ measurement of the kinetic-scale 3D anisotropy is still lacking. Moreover, to our knowledge, the investigations concerning the scale-dependency of the anisotropy, especially the scaling relations, have not been made yet. The intention of this paper is targeted to address these issues. Based on MMS measurement of magnetic field and ion velocity with unprecedented time-resolution, we applied for the first time fivepoint second-order SF to statistically quantify the 3D anisotropy of the magnetic turbulence at sub-proton scale (< 100 km). The new observational evidence we obtained, such as the empirical relations of the anisotropy scalings, can be compared with recent theoretical and numerical results, which may facilitate our understanding of kinetic scale turbulence in the space plasmas.
The paper is organized as follows. We describe the data and methods in Section 2, present an typical event with three-dimensional anisotropy in Section 3, provide the statistical results in Section 4, summarize and discusses our results in Section 5.

DATA AND METHODS
The burst mode data from four MMS spacecraft, including magnetic field (128 Hz) from FIELD instrument (Russell et al. 2016) and ion moments (6.7 Hz) from FPI instrument (Pollock et al. 2016) are used in this study. 349 MSH intervals from September 2015 to December 2019 have been selected for the statistical analysis, whereas a 10 minutes event on 4 October 2017 is presented to show a typical event with anisotropy signatures.
To quantify the anisotropy, we use the five-point second-order SF in this study. Compared with two-point SF , five-point SF is more suitable for studying spectral anisotropy at sub-ion regime (Cho 2019;Landi et al. 2019). More details of the difference between multi-point SF s are provided on APPENDIX A.
The five-point SF is the ensemble average of the squared variation ∆f from 5-point, as function of displacement l, the S The spatial variation from 5-point is measured as Where ∆f can represent perpendicular, parallel, or total magnetic field fluctuations, ... r is ensemble average over positions r. By studying the three-dimensional distribution of the S 2 (l; f ) with respect to l, the statistical shape of the eddies in turbulence can be thus inferred from the contours of S (5) 2 (l; f ). In the computation, the velocity field is used to link the time scale with space displacement as l = τ · V according to Taylor hypothesis (Taylor 1938), which has been tested in In APPENDIX B to be valid for most of our events. Here V can be adopted as the mean velocity V mean during the interval of interest or considered as the local velocity field V local , where V is interpolated on the resolution of B in the calculation. Most of previous works performed in solar wind and magnetosheath use V mean for simplicity (i.e., Chen et al. (2010Chen et al. ( , 2012; Chen & Boldyrev (2017); Wang et al. (2019b); Wu et al. (2019a)). Also, the study of electron scale magnetic field structure functions is based on V mean (Chen & Boldyrev 2017). In a recent paper by Verdini et al. (2018) the authors have considered the effects of local velocity by using V local in their analysis of velocity filed structure functions. We use five-point V local in this paper, which is defined as V local = [V (r − 2l) + 4V (r − l) + 6V (r) + 4V (r + l)) + V (r + 2l)]/16. For the small spatial scales considered here, let us rewrite l = τ · V as l = τ · (V 0 + δV ). Since τ is small, the major contribution of velocity term is from the large scale mean velocity V 0 . In other words, the resolution of V is not the key factor for interring l, since it only determines δV , whose amplitude is generally much smaller than V 0 in the magnetosheath environment under this study. To verify this point, we have performed and compared the analyses on basis of V mean and V local , and the results turn out to be nearly identical. Hence, we justify that the effects of interpolating V at the time points of B, or only considering V mean is negligible in the analysis of small-scale magnetic field structure functions.
Once the SF with respect to l is computed, it can be projected into local coordinates with respect to the local magnetic field B local as in Chen et al. (2012); Verdini et al. (2018) to study its 3D features. This coorinates allows us to compare the results with recent simulations on basis of the same reference frame (Landi et al. 2019;Cerri et al. 2019), and is consistent with previous choice of studying spectral anisotropy at various scales (Chen et al. 2012;Chen & Boldyrev 2017;Verdini et al. 2018;Wu et al. 2019a). In the Cartesian coordinates system (L ⊥ ,l ⊥ ,l ),l is along B local = [B(r − 2l) + 4B(r − l) + 6B(r) + 4B(r + l)) + B(r + 2l)]/16,L ⊥ is the "displacement direction" along δB ⊥ = B local × [δB × B local ], andl ⊥ =l ×L ⊥ . The Cartesian system can be also converted into spherical polar coordinates system as (l, θ B , φ δB ⊥ ), where θ B represent the angle between B local and l, φ δB ⊥ represent the angle between L ⊥ and the projection of l on the plane perpendicular to B local . Similar to Chen et al. (2012), angles greater than 90 • are reflected below 90 • to improve scaling measurements accuracy. Specifically, by setting the ranges of θ B and φ δB ⊥ , the SF in the three orthogonal directions can be obtained as By equating the value between pairs of S 2 (l ⊥ ), S 2 (l ), and S 2 (L ⊥ ), we could infer the anisotropy relation between l ⊥ , l , and L ⊥ .

EXAMPLE OF LOCAL 3D TURBULENCE
Here we present an example with clear signatures of 3D anisotropy. The event is observed downstream of the quasiparallel shock during 08:02:13-08:12:33 on 4 October 2017. Figure 1 shows the overview of the event. As plotted in Figure 1a-b, the magnetic field is around 6.08 (0.61, 0.06, 0.79) ± (4.1, 3.6, 3.8) nT, exhibiting numerous large directional changes while keeping its magnitude. In contrast, the ion velocity is quite stable at 288 (-0.88, 0.46, 0.07) ± (15.5, 10.2, 13.5) km/s. Figure 1c shows the instantaneous increment of the total magnetic energy δB 2 as a function of scale and time. The magnitude of δB 2 generally increase with the increase of scales and is changing intermittently with time. Similar to δB 2 , the instantaneous increments of the perpendicular energy δB 2 ⊥ and parallel energy δB 2 also exhibit the same trend with respect to spatial scale as seen in Figure 1d-e, while the former one is stronger. Figure  1f-g plots the corresponding θ B and φ δB ⊥ , respectively. Due to the rapid rotations of magnetic field, the distribution of θ B and φ δB ⊥ covers a wide range within (0, π) during the whole interval, thus allowing us to infer the 3D anisotropy of magnetic turbulence with sufficient data points. c-e) Instantaneous total, perpendicular, and parallel magnetic energy as a function of scale and time. f) Instantaneous angle between local magnetic field and space displacement vector l, as a function of scale and time. g) Instantaneous angle between L ⊥ and the projection of l on the plane perpendicular to local magnetic field. Figure 2 present the structure functions and the corresponding anisotropy scalings for the above event. The values of SF s are obtained from four MMS spacecraft, then binned and averaged. Each bin is required to have a minimum number of 200 data points to ensure reliable results as in Chen et al. (2010). For the SF of the total magnetic field energy as projected in the (l , l 2 ⊥ + L 2 ⊥ ) plane, the contours of S (5) 2 (l; B) are elongated in the parallel directions (Figure 2a), where the values at perpendicular direction are much larger than the ones in the parallel direction (i.e. S (5) 2 (l; B) at l ⊥ = 60 km is more than 100 times larger than the one at l = 60 km). This signature indicate sub-ionscale (l < 2 d i ) anisotropy with k ⊥ k , and is in agreement with Chen & Boldyrev (2017). Furthermore, we find that the contours of S 2 (l; B) are also elongated in the "displacement" direction as seen in the (l ⊥ , L ⊥ ) plane ( Figure  2b), suggesting the three-dimensional characteristics of the anisotropy. In addition to S (5) 2 (l; B), we also consider the contribution of SF by the perpendicular and parallel magnetic field, namely the S     Chen et al. (2010Chen et al. ( , 2011Chen et al. ( , 2012, which may imply the less damped state of the more anisotropic fluctuations. On the contrary, the contours of S Let us draw the attention of the reader to the point that, the sampling of SF s in the perpendicular and parallel direction is dissimilar as seen in Figure 1f, whereas the parallel data are much discretely distributed. To test whether the stationarity of the sampling will have an effect on the results of SF s, we have divided the time-series into two sub intervals and analyse the SF s separately. During Interval 1 (08: 02: 13 -08: 07: 13 UT), the measurements along parallel directions (i.e., θ B <5 • ) are less than the ones at oblique directions (θ B >50 • ). In contrast, during Interval 2 (08: 07: 13 -08: 12: 33 UT), the measurements along parallel directions are much frequent and the overall sampling are more homogeneous than Interval 1. As expected, the SF s for Interval 1 have no measurements within the range of 60 km < l < 100 km, 0 km < l 2 ⊥ + L 2 ⊥ <18 km), while the SFs for Interval 2 cover the complete wavenumber space. More importantly, the extension feature of the contours along l direction in these two sub-intervals appears quite similar as compared with the results from the whole interval. Hence, the anisotropic features of the turbulence could be viewed as stationary regardless of the interval selection. In fact, as the first step of computing SF s, we calculate the time difference between continuous sampling points rather than discontinuous sampling points. As the second step, the calculated time differences satisfying certain θ B conditions are collected together from discretely(discontinuously) distributed time points. This discontinuous collection will not significantly influence the analysis results as long as the whole time interval is statistically time stationary.
To inspect the scale-dependency of the anisotropy more precisely, we have computed the SF in the three orthogonal directions as defined by equation (3)-(5). For the 1D SF of the total magnetic field shown in Figure 2g, the relation of S 2 (l ⊥ ) > S 2 (L ⊥ ) > S 2 (l ) are satisfied at all scales as expected, thus confirming the 3D nature of anisotropy again. Moreover, this anisotropy is found to be scale-dependent. For example, at energy level of 0.01 nT 2 , the perpendicular length of the "statistical eddy, l ⊥ ∼ 6 km is smaller than the "displacement length L ⊥ ∼ 8 km, while the parallel length is much larger at l ∼ 35 km. As the energy level increase to 1 nT 2 , l ⊥ , L ⊥ , and l becomes approximately 60 km, 90 km, 150 km, respectively. The change of l ⊥ : L ⊥ : l ratio from 0.17: 0.23: 1 to 0.4: 0.6: 1 suggests that as scales increase, the anisotropy between parallel and perpendicular lengths becomes weak, while the anisotropy between two perpendicular lengths in the perpendicular plane remains almost unchanged. Setting an energy range as [0.01, 0.5] nT 2 , the SF s could be fitted by the power laws as l 1.92 ⊥ , L 1.82 ⊥ , and l 3.08 , where the standard error of the mean is 0.09, 0.08, 0.10, respectively. The power law index of the second-order structure function, g, is usually related to the power spectral index, α, by α = g + 1 (Chen et al. 2010). Hence the spectral indices in the three directions are 2.92, 2.82, and 4.08, respectively. The perpendicular spectral indicies are close to 8/3 but steeper than 7/3, which are consistent with previous findings both in the MSH and SW (Alexandrova et al. 2008;Huang et al. 2014;Matteini et al. 2017;Chen & Boldyrev 2017;Alexandrova et al. 2009). For the SF s of δB ⊥ , the trends plotted in Figure 2i are essentially the same compared with the results of δB in Figure 2g, suggesting a dominant contribution of perpendicular magnetic field fluctuations to SF s. This point agrees with the "variance anisotropy found in Chen et al. (2010) and is again supported by examining SF s of δB in Figure 2k, whose magnitudes are weaker than SF s of δB ⊥ in Figure 2g. Figure 2h displays the anisotropy relations for δB, where the blue curve represent l vs l ⊥ and the red curve represent L ⊥ vs l ⊥ . On one hand, as the perpendicular scales decrease from 1.0 d i to under 0.05 d i , the ratio of l /l ⊥ first increase and reached a maximum of ∼ 8 at ∼ 0.1 d i , whereas the anisotropy scaling obeys l ∝ l 0.67 ⊥ . Then the ratio of l /l ⊥ begins to decrease and finally approaches 1 at ∼ 0.04 d i . On the other hand, the ratio of L ⊥ /l ⊥ keeps steady around 1.5, obeying L ⊥ ∝ l 1.09 ⊥ . As expected, the anisotropy relations for δB ⊥ as shown in Figure 2j is quite similar with the relations in Figure 2h, following l ∝ l 0.59 ⊥ and L ⊥ ∝ l 1.05 ⊥ . Nevertheless, the anisotropy relations for δB as shown in Figure 2l are dissimilar, following l ∝ l 0.93 ⊥ and L ⊥ ∝ l 1.01 ⊥ . Figure 2. The SF s in the 2D plane, together with the 1D SF s and their anisotropy scalings. The left two columns include: (1) 2D SF s as a function of (l , l 2 ⊥ + L 2 ⊥ ), and (2)SF s as a function of (l ⊥ , L ⊥ ). The right two columns plot respectively: (1) 1D SF as a function of l ⊥ , L ⊥ , and l , and the relations between l , l ⊥ , and L ⊥ . For each panel, the first, second, and third rows represent SF s of the total, perpendicular, and parallel magnetic field, respectively.

STATISTICAL ANALYSIS OF THE ANISOTROPY SCALINGS
In this section, the sub-ion-scale anisotropy relations are investigated comprehensively based on a statistically survey of 349 intervals during 2015-2019, when MMS instruments was in burst mode. These intervals are tagged as "magneotosheath on the MMS science data center website. In addition, to avoid the influence of shock or magnetopuase, the ion and electron energy spectrogram have been checked by eye to make sure they exhibit typical broadband MSH signatures and are time stationary. As a result, 349 intervals with an average duration of ∼ 5.8 minutes have been selected. Figure 3 shows the histograms for the events duration and plasma parameters, together with their mean value and standard deviations. The events duration, proton beta β p , ion inertial length d i , and proton gyroradius ρ p , cover the range of [60, 1680] s, [0.3, 80], [15, 130] km, and [45, 370] km respectively. The mean value of temperature anisotropy T p⊥ /T p is 1.05. The distribution of mirror mode threshold Σ mirror = T p⊥ /T p − 1/β p⊥ − 1, is also shown in Figure 3e. Since most of the values are negative, the influence of mirror instability is not strong in our database.  Regarding the anisotropy of δB and δB ⊥ in the perpendicular plane, Figure 4d-e show that for most of the data, although L ⊥ > l ⊥ , the anisotropy level is stable since the scalings are close to 1. In addition, the empirical relation L ⊥ = b 0 · l β ⊥ for the median value are fitted as L ⊥ = 1.58 · l 1.08±0.01 ⊥ , L ⊥ = 2.03 · l 1.13±0.02 ⊥ at [0.02, 1.0] d i , respectively. Lastly, the anisotropy for δB vanishes and follows a nearly isotropic relation of L ⊥ = 1.15 · l 1.01±0.003 ⊥ .

CONCLUSION AND DISCUSSION
In this paper, we have conducted a statistical survey of the sub-ion-scale anisotropy of the turbulence in the Earths magnetosheath. By measuring the five-point second-order SF s of the magnetic field, the three-dimensional structures of the turbulence have been quantitatively characterized. Specifically, the three characteristic lengths of the eddies are found to roughly satisfy l > L ⊥ > l ⊥ in the local reference frame defined by Chen et al. (2012). As for the scaledependency of the anisotropy inferred from SF s of the total magnetic field, (1) the parallel−perpendicular anisotropy as revealed by the ratio of l /l ⊥ , shows an increase trend towards small scales and obeys a scaling of l 2.44 · l 0.71 ⊥ between 0.1 d i and 1 d i , then it decreases and obeys l 170 · l 2.45 ⊥ between 0.04 d i and 0.08 d i .
(2) the anisotropy in the perpendicular plane as revealed by the ratio of L ⊥ /l ⊥ , is generally weaker than the ratio of l /l ⊥ . Moreover, this anisotropy is scale-invariant, displaying a scaling of L ⊥ 1.58 · l 1.08 ⊥ . Interestingly, the parallel−perpendicular anisotropy tends to become increasingly isotropic when approaching both large scales ∼ 4 d i and small scales ∼ 0.04 d i (Figure 4 a-c). This large-scale isotropy may reflect similar structures as in the isotropic solar wind reported recently (i.e., Wang et al. (2019b); Wu et al. (2019a,b)), while the small-scale isotropy has not been reported before to our knowledge. Possible explanations for such isotropy include the weakening of perpendicular cascade and the influence of ion cyclotron waves. Indeed, there are a few events with l ⊥ > l in the database, where the existence of ICW has been confirmed by checking the polarisation state of the fluctuations. In a few other events, we also find coexistence of ICW and 2D l > l ⊥ structures through inspecting the SF s in the 2D (l , l ⊥ ) plane.
The scaling of l ∝ l 0.71 ⊥ observed in this work is different from traditional KAW theory of l ∝ l 1/3 ⊥ , but is close to the theoretical prediction of l ∝ l 2/3 ⊥ from Boldyrev & Perez (2012), Boldyrev & Loureiro (2019) and simulation from Cerri et al. (2019). It also corresponds to α = 1.13 in the framework of the model by Landi et al. (2019). Hence a modified CB premise may be needed to understand the kinetic turbulence in magnetosheath. We note that, for most of the events, the week anisotropy in the perpendicular plane is inconsistent with results of (Boldyrev & Figure 4. Statistical analysis of the 3D anisotropy scalings. The first, second, and third column shows respectively the anisotropy of δB, δB ⊥ , and δB . The light grey curves represent the superposition of the statistical results, the median value and fitted results are overplotted in bold solid and dotted lines. The dashed lines represent specifically, the reference scalings with slope of 1/3, 2/3, 1 at top panels, and 4/2, 3/2, 1 at bottom panels. The histogram of the anisotropy scalings at 0.04 di -0.08 di (light grey), and 0.1 di -1.0 di (dark grey) are inserted in top panels. Similarly, the histograms of the anisotropy at 0.02 di -1.0 di (light grey) are inserted in bottom panels.
Loureiro 2019), which predicts a much stronger anisotropy and a steeper scaling of L ⊥ ∝ l 3/2 ⊥ or L ⊥ ∝ l 2 ⊥ , depending on different current sheet configurations for the tearing instability. Capturing the active signatures of current sheet disruption/reconnection (i.e., (Mallet et al. 2017;Loureiro & Boldyrev 2017)) from in-situ observation is a challenging task, but will contribute to understand its effects on the anisotropy.
The spectral anisotropy of kinetic plasma turbulence is believed to be associated with dispersion and intermittency effects (Zhao et al. 2016;Landi et al. 2019). For example, the anisotropic scalings are different below and above ion cyclotron frequency and also differs for sheet and tube like turbulence (Zhao et al. 2016). To illustrate the possible connection between intermittency and spectral anisotropy, we specifically compare the results from two events. Figure  5 shows the excess Kurtosis and anisotropy relation for the total magnetic field. The excess Kurtosis is defined as K(l) = S 4 (l)/S 2 (l) 2 − 3, where S 4 (l) is the fourth-order structure function. In both panels of Figure 5, the solid lines represent the results from event 1, which is recorded on 4 Oct 2017 and is used as our example of spectral anisotropy in section 3, while the dash-dot lines represent the results from event 2, which is recorded on 24 Dec 2017. As plotted in three directions l ⊥ , L ⊥ , l , the value of K(l) is around zero at large scales, meaning the roughly Gaussian distribution of the fluctuations. Toward small scales, K(l) displays an increase tendency before it drops down. Such non-Gaussian statistics (K(l) > 0) confirms the presence of intermittency in the magnetosheath, while the scale-dependent profile of Kurtosis is similar with solar wind observations (He et al. 2019). For the kinetic scale parallel-perpendicular anisotropy, we find that it can be considerably affected by the intermittency. As shown in Figure 5, the stronger the intermittency (see the larger Kurtosis of the solid lines in the left panel), the stronger the anisotropy level (see the larger vertical distances between the solid blue curve and the grey line in the right panel). This phenomenon is consistent with previously observations at large scales, which emphasize the key role of intermittency in generating the spectral anisotropy (i.e., (Wang et al. 2014;Pei et al. 2016;Yang et al. 2018)). For the anisotropy in the perpendicular plane, the anisotropy levels (see the two red curves in the right panel) of these two events are much smaller as compared with the parallel-perpendicular anisotropies. The reason for such weak anisotropy still remains to be explored. We note that it has been proposed that the axisymmetric, 2D (k ⊥ > k fluctuations can be observed with non-axisymmetric features in the spacecraft frame due to a sampling effect (i.e., Alexandrova et al. (2008); Turner et al. (2011);Lacombe et al. (2017) and Matteini et al. 2020, submitted). Therefore, the spectral anisotropy (non-axisymmetric) in the perpendicular plane needs be cautiously interpreted with such effects being quantitively explored in the future. Lastly, we note that both of the events still have non-Gaussian fluctuations (as seen in the non-zero values of the Kurtosis), hence the absence of isotropic l ≈ L ⊥ ≈ l ⊥ relation is not contradictory to previous studies, which found isotropy when the intermittency is removed (i.e., (Pei et al. 2016)). We plan to conduct a much comprehensive analysis to understand how intermittency influence spectral anisotropy in a future work, particularly focusing on comparing the role of different coherent structures on the anisotropy (i.e., 2D tube-like vortices in Wang et al. (2019a) or 1D current sheets. The differences between multi-point second-order structure functions of the total magnetic field are compared here, where the two-point and three-point structure functions are defined as S 2 (l; f ) = |f (r+l)−f (r)| 2 r , and S 2 (l; f ) = 1/3 |f (r − l) − 2f (r) + f (r + l)| 2 r , respectively. As seen in the left panel of Figure 6, the trend of three-point and fivepoint SF tend to agree with each other, whereas the slope of two-point SF are relatively flatter, especially towards small scales. We also compute the "equivalent spectrum defined as S 2 (k)/k and compare the results with power spectral density (PSD) from Fourier transform. Again, it is found that within [0.05, 1] km −1 of the right panel, the slope of the spectrum based on five-point SF is around -2.86, which is similar with the three-point result of -2.80 and the slope of PSD around -2.94, while the slope based on two-point SF is only -2.55. Hence it is demonstrated that in order to capture the spectral characteristics of the turbulence at sub-ion regime, the use of multi-point (>2) structure functions are preferred.  (2019)). The validity of Taylor hypothesis for all the events is checked by comparing the structure function of magnetic fluctuations in two ways (Chen & Boldyrev 2017): one is to calculate the structure function from single-spacecraft measurements by assuming Taylor hypothesis, and the other is to calculate the structure function based on direct spatial differences between measurements from six pairs of MMS spacecrafts, which are separated by certain inter-distances between them. Figure 7 shows the statistical results of the first-order structure function as a function of scale. As represented by different colour for each individual event, the results based on Taylor hypothesis (curves) are close to the results from direct spatial measurements (crosses) at 5 km < l < 200 km. Therefore, the use of Taylor hypothesis in our analysis have been proven to be reasonable. We note that our results at small scales are in agreement with recent demonstration of Taylor hypothesis being valid down to kd e ∼ 1 (Chen & Boldyrev 2017). In addition, the Taylor condition are known to be better satisfied at relatively large perpendicular wavenumbers especially when fluctuations are sampled along the perpendicular direction (Chen & Boldyrev (2017) and references therein), thus the presence of spectral anisotropy (k ⊥ > k ) in our events is also in favour of the Taylor assumptions. Figure 7. Test of validity of the Taylor hypothesis. The vertical axis represents the normalized magnetic fluctuation amplitudes obtained from first-order structure function. The curves are from the time differences of measurements from individual spacecraft under Taylor assumption, while the crosses are based on the spatial differences between measurements from the six pairs of MMS spacecraft.