Anisotropy of the Spectral Index in Ion Scale Compressible Turbulence: MMS Observations in the Magnetosheath

Turbulence in the Earth’s magnetosheath at ion kinetic scales is investigated with the Magnetospheric MultiScale (MMS) spacecraft. The multi-point measurements allow the three dimensional power spectra in wave-vector space to be determined. Previously the three dimensional structure of fluctuations in the magnetic field and density (using spacecraft potential as a proxy) were possible with Cluster. However, using the excellent time resolution data set provided from both the Fluxgate Magnetometer (FGM) and the Fast Plasma Investigation (FPI) on MMS the spectra can be determined for a number of different parameters such as ion velocity, and ion temperatures parallel and perpencidular to the mean magnetic field directions. The spectra for different fluctuations show similar features to one another such as a strong power anisotropy with respect to the mean magnetic field direction, such that the energy decays faster in the direction parallel to the mean magnetic field than the perpendicular direction. A weak non-gyrotropy is also seen in the direction of the bulk velocity similar to what has been seen in magnetic field fluctuations with Cluster at ion kinetic scales in the solar wind. Velocity fluctuations are shown to be the most anisotropic. The density and temperature fluctuations exhibit similar anisotropies but are much weaker in comparison.

Turbulence in the Earth's magnetosheath at ion kinetic scales is investigated with the Magnetospheric MultiScale (MMS) spacecraft. The multi-point measurements allow the three dimensional power spectra in wave-vector space to be determined. Previously the three dimensional structure of fluctuations in the magnetic field and density (using spacecraft potential as a proxy) were possible with Cluster. However, using the excellent time resolution data set provided from both the Fluxgate Magnetometer (FGM) and the Fast Plasma Investigation (FPI) on MMS the spectra can be determined for a number of different parameters such as ion velocity, and ion temperatures parallel and perpendicular to the mean magnetic field directions. The spectra for different fluctuations show similar features to one another such as a strong power anisotropy with respect to the mean magnetic field direction, such that the energy decays faster in the direction parallel to the mean magnetic field than the perpendicular direction. A weak non-gyrotropy is also seen in the direction of the bulk velocity similar to what has been seen in magnetic field fluctuations with Cluster at ion kinetic scales in the solar wind. Velocity fluctuations are shown to be the most anisotropic. The density and temperature fluctuations exhibit similar anisotropies but are much weaker in comparison.

INTRODUCTION
Turbulence is a phenomenon characterized by disordered fluctuations in several of the fluid's parameters over a large range of time and length scales. For a neutral fluid this might include density, velocity, and temperature, however in a space plasma there are also fluctuations in electromagnetic fields due to a very high conductivity [1,2]. Typically most of the research on the topic of in situ plasma turbulence have been performed using magnetic field data as the data are often operationally simpler to obtain with high time resolution. To obtain a more complete understanding of the turbulent fluctuations, measurements are required for parameters other than the electromagnetic fields. While density can be obtained from spacecraft potential (e.g., [3][4][5]), other plasma measurements such as velocity and temperature require a direct plasma measurement. Typically plasma instruments are mounted looking in one direction on a spinning spacecraft and use the spin to obtain data azimuthally. Thus, the time resolution is limited to the spacecraft spin, which is typically not fast enough to resolve ion kinetic scales.
The majority of plasma turbulence studies have employed single spacecraft measurements, where spatial information is obtained by assuming Taylor's hypothesis [24], where the turbulent fluctuation is a assumed to vary slowly with respect to the measurement time. By analyzing intervals with different orientations of the magnetic field with respect to the bulk flow direction, the three dimensional structure is inferred. Using this approach single spacecraft observations have revealed correlations between solar wind measurements in the directions perpendicular to the magnetic field are longer than in the parallel direction [25] giving a "Maltese Cross" pattern showing the dominance of perpendicular wave-vectors. Furthermore, when the data are classified as fast (typically above 500 km/s) the opposite is true [26] with parallel wave-vectors dominating. This same pattern has been observed in variables other than magnetic field by Smith et al. [27], where velocity, temperature, and density show similar correlations as the magnetic field. Recently Wang et al. [28], used the self correlation technique and revealed with intervals to reveal similar anisotropies to Dasso et al. [26] for intervals of length 1-2 days, but become more isotropic as the intervals become smaller i.e., when looking at smaller scales in the inertial range.
The power of magnetic fluctuations in the solar wind at large inertial scales [where a magnetohydrodynamic (MHD) description is valid] have been found to be generally smaller along the magnetic field direction with respect to the perpendicular direction and the spectrum in this direction is steeper with a spectral index of −2 [9,29] compared with an index of −5/3 in the perpendicular direction. Structure function analysis has also been performed on magnetic field data revealing anisotropic power that also evolves toward smaller scales [22], and that fluctuations in the magnitude of the magnetic field (a proxy for the compressible fluctuations) are more anisotropic than the trace fluctuations in the fast solar wind.
At smaller scales often plasma data are usually not available as the instruments lack the necessary time resolution. This limitation allows the study of plasma fluctuations only in the low frequency inertial range where a fluid description is still valid. At higher spacecraft frame frequencies (denoted by subscript sc) above around f sc 1 Hz the fluctuations become comparable to proton gyration frequencies. Above these frequencies the protons cannot follow the magnetic field any longer and become demagnetized, while electrons remain magnetized and can still follow the magnetic field due to their smaller gyroradius. This region is often marked by a break in the power spectral density and a steepening of the spectrum [5,[30][31][32][33][34][35]. The location of the break has been measured to be fairly independent of plasma β [33,34] and varies with heliocentirc distance [32]. Both of these observations can be explained by the break corresponding to the scale of the ion cyclotron resonance. Additionally the location of the break has been observed by single spacecraft observations to be independent of the angle between the magnetic field and the bulk velocity [35]. This was interpreted by Duan et al. [35] due to the ion diffusion region (where ions decopule from electrons) being approximately isotropic in wave-vector space. The scales smaller than the spectral break is often termed the dissipation/dispersion range, or the ion kinetic range, and a fluid description is no longer valid.
In some conditions Taylor's hypothesis may break down, should turbulent fluctuations become very dispersive, when bulk speeds are low, or different modes appear at once in the plasma [36][37][38]. The magnetosheath is an especially interesting plasma as it typically has a higher magnetic field strength, a lower bulk velocity, and has a much larger compressibility than the solar wind. The lower speed and larger fluctuation amplitudes make the breakdown of Taylor's hypothesis in this region more likely when juxtaposed with the super Alfvénic solar wind. To overcome the limitations imposed by Taylor's hypothesis two other approaches to understand the structure of turbulence are possible which are through either multi-point measurements, or from direct numerical simulation.
Different numerical simulation schemes have been used to investigate the structure of turbulence, however recently an expanding box MHD approach can also include expansion effects present in the solar wind allowing the large scale three dimensional structure to be simulated [19,39]. The other approach is to use multi-point measurements, which are possible with Cluster and MMS. These techniques rely on differences or correlations between spacecraft pairs e.g., multi spacecraft structure functions/cross correlations [40][41][42] where gyrotropy is assumed. These techniques have also revealed anisotropic power and spectral indices parallel and perpendicular to the mean field direction. Another approach is to use multi-point signal resonator technique [13] which assumes a plane wave geometry of the fluctuations and make use of the phase delay between measurement points [5,13,43] but does not make an assumption of gyrotropy. However, these studies performed with Cluster have mostly been performed using magnetic field data, or spacecraft potential [5,18]. This study will expand on the previous work done by Cluster in the inertial and kinetic ranges, and use the Magnetospheric MultiScale mission [44] and its exceptionally high time resolution plasma data from the Fast Plasma Investigation [45] (FPI) instrument to characterize the structure of turbulent fluctuations in the transition from fluid to kinetic scales in the Earth's magnetosheath.

DATA/METHODOLOGY
We use data from the MMS spacecraft [44]   of studying magnetic reconnection. This is ideal for investigating the scales near proton gyration/inertial scales. The spacecraft were in the dusk side flank of the magnetosheath downstream of the quasi-perpendicular shock. A summary of the interval is shown in Figure 1 which shows the magnetic field, ion velocity, ion density and ion and electron temperatures and mean parameters of the interval are given in Table 1.
The multi-point signal resonator technique [13] (MSR) will be used to analyze the different measurements given in Figure 1. Whereas, previously this technique (or it's predecessor k-filtering/wave telescope) has been applied to electromagnetic [46], magnetic fields [14], and to density derived from spacecraft potential [5]. In this study we will investigate magnetic, velocity, temperature, and density fluctuations. For direct comparisons to the plasma data the magnetic field data will be resampled to the ion measurement time tags. The MSR technique relies on weak time stationarity, and spatial homogeneity of the signal. The signal seems fairly homogeneous in terms of the mean value throughout the interval, however after 14:03:30 the fluctuations seem a little less Alfveńic. This however makes up only a minority of the overall signal and is unlikely to have a large effect on the results based on testing with slightly different intervals than shown here. It should also be noted that a fine balance needs to be struck between the need sufficient data points for ensemble averaging and the weak stationarity of the signal. Especially when intervals such as this one where the MMS separations are large enough to study the ion kinetic range are extremely rare. As an additional test we applied the same analysis to the full resolution magnetic field (where more data points are available than the ion data) of the shorter interval up to 14:03:30 and the method yielded similar results for the magnetic field. Thus, we can be confident that the analysis is justified for this interval. The MSR technique also supposes that the signal can be described mathematically as a superposition of plane waves with random phases in the spirit of Fourier analysis, and a small component of incoherent noise. Essentially the technique uses the time series sampled at each spacecraft and the signal can be filtered for a frequency ω sc using a Fourier transform, and for wave-vector k using the multiple measurement points. Thus, only power related to a plane wave with frequency ω sc and wave-vector k is transmitted through the filter, thus by investigating a number of wave-vectors a distribution of power in wave-vector space can be estimated at a given frequency P(ω sc , k). It is important to note that this does not require that the different plane waves that the signal can be decomposed into correspond to any particular linear wave solutions of the Vlasov equation.
In the case of using magnetic or velocity data, three components at each spacecraft can be used as an input giving a total of 12 time series, however for magnetic field data the filter can be further refined by enforcing the solution to conform to the divergence free condition for magnetic field which is termed a constraining matrix [47]. The application of the method to density is detailed in Roberts et al. [5], where there are only four time series input and similarly to the velocity case no constraining condition can be imposed. The unique capabilities of MMS allow it to be used on velocity, temperature fluctuations as well as a direct measurement of the density. In this study we will use the DC magnetic field from the fluxgate magnetometers (FGM) [48] which have a sampling rate of 128 Hz in burst mode and the ion plasma measurements from the Fast Plasma Investigation's Dual Ion Spectrometers [45] (FPI-DIS) which has a rate of 6.6Hz. As previously mentioned the magnetic field data is resampled to the ion time tags.
The method is subject to some limitations; the smallest scale that can be investigated is limited by the mean inter-spacecraft distance k max = π/ d which is the primary driver for the choice of interval. The value k max defines a cube in wave-vector space that extends from −k max to +k max , such that the length of one of the sides of the cube is 2π/ d . The k max length scale is related to a timescale for an advected structure giving an upper frequency bound of f max = V i /2 d . Where we take the bulk speed to be the ion bulk speed. Conversely the large scale limit is set to k max /25 when the error of determining a wave-vector becomes larger than 10% for a simulated plane wave [5,49]. The technique also assumes weak stationarity and that the fluctuations can be described as a superposition of plane waves and incoherent noise. The method has been tested for a signal made of intermittent bursts of activity [50], where it is shown that coherence in the signal does not affect the ability to resolve incoherent features. Furthermore, spatially repeating coherent structures (with a Poisson distributed size and spacing) can be recovered, and resemble wave-packets that have a random phase. Additionally a test was performed using both the MSR method and multi-spacecraft timing, which are based on different assumptions yielded similar results for an interval where several different intermittent coherent structures were identified [51].
In this work we will analyze the results from the MSR method in two different ways; Firstly P(k) will be obtained by integrating P(ω sc , k) with respect to the plasma frame frequency to investigate the anisotropies in the power distribution. This integration is performed between the limits f min = 0.06, f max = 1 Hz which come from the spacecraft mean separation ∼ 140 km and the bulk flow speed. This will give a measurement in the power in wave-vector space, where we will for each measurement quantify the possible anisotropies and agyrotropies present. Secondly, we will also reduce the three dimensional spectra to spectra along one direction to investigate how the spectral index in the dissipation/dispersion range varies with the angle from the magnetic field.

One Dimensional Spectra
We begin this section by investigating the typical one dimensional analysis usually performed by investigating the Fourier spectra of the various different plasma parameters. Different Fourier spectra are shown in Figure 2, for magnetic field, ion velocity density and parallel and perpendicular ion thermal speeds. The frequencies corresponding to the proton gyroradius and inertial lengths f 87 Hz are indicated in orange assuming a mean bulk flow indicated in Table 1. The velocity and magnetic spectra are fitted with two power laws, the density and perpendicular thermal speed spectra are fitted with three power laws while the parallel thermal speed is fitted with one power law. The error on the spectral indices is obtained from the residuals of linear least squares fitting for log power against log wave-number. The break frequency is found by fitting the two power laws from opposite sides of the spectral break and then finding the intersection of the two lines. This procedure is done twice for the density/perpendicular thermal speed spectra as there are two break locations. Figures 2A,B show the spectra of the magnetic field (in Alfvén units) resampled to the ion velocity time tags and the ion velocity respectively. At large scales a spectral index close to -5/3 is obtained for the magnetic field with the velocity spectra being noticeably shallower, closer to -3/2 as is often observed in the solar wind and in the magnetosheath [52][53][54]. The inertial range is not always seen in the magnetosheath especially in the outer magnetosheath as the interaction of the solar wind with the Earth's bow shock destroys the correlations in the inertial range and results in a 1/f spectrum which transitions straight to the ion kinetic range without having time for an inertial range to develop.
In this case the observations are taken in the inner magnetosheath such that a well-developed inertial range is seen and is followed by a break and then a steepening with the ion velocity being steeper than the magnetic field [55,56]. The magnetic field fluctuations dominate the velocity fluctuations in this case which may be due to the development of 2D structures such as the Alfvén vortex at inertial scales [51,[57][58][59][60] or current sheet or flux rope like structures [61][62][63]. At a spacecraft frame frequency near ∼ 0.2 − 0.4 Hz there is a spectral break which is far from both the Taylor shifted inertial and Larmor scales. In the solar wind for extreme values of β the break location (e.g., [64]) has been shown to agree better with the larger of the two scales, although for typical values of β in the solar wind there is not much effect [34]. Alternatively the break in the solar wind shows good agreement for a variety of heliocentric distances with the scale expected for cyclotron resonance [65] Hz which we denote in Figure 2 as a purple vertical line. For the case presented it seems that the measured break point is more closely related to that expected for cyclotron resonance. Interestingly, for the density spectrum the first break agrees well with the f c and the second break agrees well with the other break scales. A small bump is seen in this range for the perpendicular thermal speed, while no clear break can be found in the parallel thermal speed and only a single power law is fitted. The bump in the perpendicular thermal speed spectra could be due to cyclotron resonance (e.g., [33,66,67]) which could act at this range of scales and the wave particle interaction would be expected to heat protons in the perpendicular direction. However, no signature is seen in the trace magnetic spectra but a flattening is seen in the density spectra, which could suggest that a small scale compressible process is active. The time series in Figure 1 show that there is an anti-correlation between density and temperature suggesting that the compressible fluctuations exhibit pressure balance. One possibility is that slow waves could be responsible (e.g., [68,69]). Compressible slow waves are heavily damped in a plasma such as the magnetosheath due to the very high ion to electron temperature ratio. However, they are damped proportionally to k thus could exist as pressure balanced structures k = 0 or highly oblique kinetic slow waves k ⊥ ≫ k . There is some evidence from numerical simulations of the magnetosheath behind the quasi parallel shock that suggests that slow waves can exist in these conditions [69].
At kinetic scales the dominance of the magnetic spectra could be due to kinetic Alfvén wave like fluctuations [56,70], rather than kinetic slow waves. This suggests that should slow waves be responsible for the bump in the perpendicular temperature and the density they are not dominant at smaller scales. The gray trace show the estimation of the FPI velocity noise floor [55], showing that noise becomes significant near f sc = 1.5 Hz in the ion data. Figure 2C shows spectra of the ion density, contrary to the magnetic field and velocity spectra a single clear break is not identified. Rather there exists a transition region which begins near the expectation for cyclotron resonance where the spectra flatten before it steepens near the inertial/Larmor scales. The transition in the density scale has been often observed in the solar wind and has been modeled as being due to the presence of compressible slow waves at large scales and kinetic Alfvén waves at small scales [71,72]. The location where these two phenomena exist is dependant on the plasma beta with a smaller beta giving a larger transition region. However, it is not clear whether this is applicable for the magnetosheath, as other phenomena such as mirror modes are much more common in the magnetosheath and contribute significantly to the compressible power at inertial scales.
The nature of the thermal speed spectra of plasma turbulence is the subject of recent debate [73]. With high time resolution measurements of the solar wind only available from the Faraday Cup on Spektr-R. These spectra showed similarities to the spectra of the velocity fluctuations [74,75] although it is argued by Gogoberidze et al. [73] that effects due to anisotropic temperatures and the measurement from the Faraday cup are misleading. In their work they argue that the proton thermal velocity should have a shape more similar to the compressible fluctuations, or the trace magnetic fluctuations. The thermal speed spectra are shown in Figures 2D,E. The perpendicular thermal speed spectra somewhat resembles the compressible fluctuations and in the same region as a flattening is seen in the density spectrum a small bump is seen at the same range.
While we have presented the Fourier power spectra of several different parameters here the key limitation is that these are along a single path of the spacecraft through the plasma (Taylor's hypothesis). But power distributions and spectral indices are anisotropic and it is possible that spectra of different parameters could resemble the spectra taken along a certain direction. We will now use the multi-spacecraft capabilities of MMS to explore the different spectra in three dimensions.

Three Dimensional Spectra
In this section we move on to multi-spacecraft analysis of the fluctuations using the MSR technique. Turbulence spectra show the power is concentrated at low frequencies (and wavenumbers). However to get a full distribution of the power in wave-vector space it is insufficient to only consider the low frequency Fourier modes as their power corresponds to power at small wavenumbers. Therefore, an integration is needed in frequency to contribute the powers at higher wave-numbers. Thus, after obtaining a four dimensional power P(ω sc , k) at each spacecraft frequency we integrate between spacecraft frame frequencies of 0.06-1Hz to obtain the power distribution in wavevector space P(k). One of the advantages of this technique is that Taylor's hypothesis is not invoked, so there are no concerns about the validity of Taylor's hypothesis when integrating higher frequencies. Figure 3 shows the three dimensional power distributions obtained through the application of the MSR technique to various different measured parameters. These distributions have been integrated over the third direction to give a two dimensional representation of the vector quantities, the magnetic field (a), the velocity (b). Meanwhile Figure 4 shows the scalar quantities density (a) and parallel and perpendicular ion temperatures (b,c). The co-ordinate system is the mean field aligned system with the e = B 0 /|B 0 |, e ⊥1 = e × V i /|V i |, and e ⊥2 = e × e ⊥1 such that the bulk flow is primarily in the −e ⊥,2 direction.
All spectra show similar shapes to one another but there are slight differences. To quantify the relative levels of anisotropy and agyrotropy in the power distributions of different powers we use eigenvalues of the stress tensor where the i,j elements of the tensor are defined in Equation (1).
Where the angled brackets denote the centre of the power distribution. This metric has often been used to quantify the geometrical properties of numerically simulated velocity distribution functions (e.g., [76]). The eigenvectors of this stress tensor give three orthogonal components [e 1 , e 2 , e 3 ] which are aligned with the major axis and two perpendicular axes where the corresponding eigenvalues satisfy λ 1 > λ 2 > λ 3 . The eigenvalues can be used to define an anisotropy index λ 1 /λ 3 and a non-gyrotropy index λ 1 /λ 2 . For a sphere all eigenvalues would be equal and for a gyrotopic cigar like shape which is elongated along one axis the two smaller eigenvalues would be equal. The results in Table 2 demonstrate that not only that magnetic fluctuations are anisotropic but all parameters investigated are anisotropic and show similar distributions. We note that the power distributions of the vector fields are more irregular than those of scalar quantities, this may be due to the number of components put into the MSR method i.e., 12 vs. 4.
The most anisotropic component is the velocity. This perhaps reflects the differences in the one dimensional spectra as proton kinetic effects become important and the velocity spectrum steepens significantly compared to the magnetic spectra [56]. Meanwhile the spectra of other parameters are much less steep than the ion velocity. The scalar components are less anisotropic than the both magnetic and velocity fields which is different to cases in the solar wind where the compressible component can be more anisotropic [22,56]. In the fast solar wind at inertial scales this has been interpreted to be due to slow waves which are damped proportionally to k , thus only the fluctuations with the most perpendicular wave-vectors can survive, while the incompressible component contains Alfvén waves which do not have such a sensitivity to the propagation angle. However, in the magnetosheath the compressible component is unlikely to be due to kinetic slow waves as they are more severely damped due to the plasma conditions i.e., high ion to electron temperature ratio (e.g., [77]). The dominance of the magnetic field spectra over the velocity spectra suggests that for this interval slow waves are not dominant at kinetic scales, but could explain the bump in the perpendicular temperature and the flattening in the density spectra. This could be due to damping before they can cascade to the kinetic range. However, simulations do suggest that slow waves can exist in this environment at kinetic scales (e.g., [69]). One possibility is that pressure balanced structures could exist at kinetic scales as they are undamped. This interpretation has a weakness, as the compressible three dimensional spectra are less anisotropic than the magnetic and velocity components whereas pressure balanced structures are expected to be highly anisotropic. We suggest that something else must contribute to the compressibility, perhaps structures such as mirror modes which are quasi aligned with the magnetic field giving a smaller anisotropy than the magnetic components.
The similarity in the shapes of the distributions echoes previous work at large scales by Smith et al. [27], where similar shapes in the correlation lengths were found when using magnetic, velocity, density, and temperature. There are also suggestions here of non-gyrotropy which seems to be larger when the anisotropy is larger. The reason for this is not fully understood, some hypotheses include sampling effects due to the velocity direction [78], a remnant of a large scale effect due to the bow shock, or a preferred cascade direction. A comparison with a simulation with a known distribution of power in wave-vector space should be performed to confirm or refute the sampling effect described by Turner et al. [78]. To be able to fully understand the causes of this anisotropy more intervals should be analyzed in the magnetosheath an the solar wind and comparisons made with numerical simulations which incorporate other effects which can cause second preferred directions to appear. There is also a weak asymmetry with respect to the e direction showing more power k < 0.

Anisotropy of Spectral Index
To further investigate the anisotropy we reduce the full distributions to a one dimensional spectrum. We focus only on the spectral index of the dissipation scale as inertial scales have too few points in wave-vector space for an accurate determination of the spectral index due to the interspacecraft distances being too small. The distributions are converted to cylindrical coordinates and the perpendicular direction is defined as the integration over the azimuthal direction, the distribution is then folded to give a spectrum in (k , k ⊥ ). Finally this is reduced to a one dimensional spectrum by taking a 1D cut through the distribution.
One dimensional spectra in directions parallel and perpendicular to the mean magnetic fields in Figures 5, 6, and in the flow direction in Figure 7. Further details of the normalization of both the MSR spectra and the conversion and normalization of the spectra in Figure 2 are given in the Appendix. It is important to note that an incorrect normalization was displayed in Roberts et al. [17] where the spectrum was divided by k rather than k and the total power was missing. Figure 5 shows the one dimensional cut along the parallel direction, for the five different quantities in this study. The scales FIGURE 4 | Two dimensional representations of the three dimensional power distribution P(k) for scalar quantities density (a), ion parallel thermal speed (b), and ion perpendicular thermal speed (c). The anisotropy is quantified as the largest eigenvalue over the smallest value. The nongyrotropy is quantified as the maximum over the intermediate eigenvalue. The orientation of the power distribution with respect to the mean magnetic field direction is calculated using the largest eigenvector and expressed in degrees θ Max = cos −1 e 1 · e .
corresponding to the Larmor radius and inertial length are shown in orange and the scale corresponding to the cyclotron frequency is shown in purple similarly to Figure 2. Figure 5 shows an interesting feature that there is an enhancement in the parallel power and in the perpendicular thermal speed. This could be an indication of ion cyclotron waves, which cannot be seen in the time series as the streamwise wave-vector is far from parallel. This signature is also at the low wave-number end which may be less accurate [49]. Confirmation could be made by investigating another interval with larger separations so that the scales can be studied with more accuracy. Future work should seek to perform comparisons with time intervals where the magnetic field is along the velocity direction where the signature would be able to be seen in the sapcecraft frame Fourier spectra. There are several results in the fast solar wind that show signatures of ion cyclotron waves for intervals with radial fields. However, it is not clear whether cyclotron waves are present but can only be seen for times when the magnetic field is pointing radially [10,79]. Figure 6 shows the same plots but for the perpendicular direction, while Figure 7 shows the reduced one dimensional spectrum for an angle of 50 • which is approximately along the bulk velocity direction. This is plotted with the spacecraft frame spectral density (corresponding to the PSD in Figure 2) with the wave-number determined by assuming Taylor's hypothesis k = 2πf /V i . Both vector spectra show reasonable agreement with the best agreement with the velocity. The results here are different to Figure 2E in [17] as in that paper the comparison with the flow direction spectra was made by averaging in the Cartesian coordinates for comparing the two perpendicular components rather than a one dimensional cut. Therefore, it is likely that the power there did not match as well due to the stronger agyrotropy which was averaged over. Additionally here a 1D cut is used through the spectra. At smaller scales the spectra differ possibly due to the Taylor's hypothesis losing validity for compressible components or due to spatial aliased power contributing at smaller scales. Spatial aliasing might be more prevalent for the scalar quantities when compared to the spectra of the time series as there are fewer inputs into the method. Aliasing is also suppressed by the constraining condition of the magnetic field. This explains why the power from the MSR method is much higher for the scalar quantities at small scales. To have a definitive conclusion about the application of the method to scalar components and to understand the agyrotropy comparisons of the method with numerically simulated data will be performed in the future. The results from Figures 5-7 are summarized succinctly in Figure 8, where the spectral indices are plotted as a function of angle from the magnetic field. The spectra are fitted at several angles from the mean magnetic field direction (diamonds) and are compared to the value for the one dimensional spectra (circles) in Figure 2. The two vector spectra agree well with the spectral indices measured in frequency space apart from the density and perpendicular thermal speed most likely due to the nature of the spectra having a transition region/bump, or due to the lack of a constraining condition/fewer time series inputs as mentioned previously. The comparison between the spectral indices in frequency space and wave-number is difficult to make in some cases. This is because the spectral break location varies between each parameter as observed in Figure 2. However, the wave-vector range available is controlled by the spacecraft separations which remains the same. Therefore, in some cases there may be more Fourier modes present in the inertial range (or kinetic range) for different parameters due to the different location of the break. This is seen in the density spectrum and the perpendicular thermal speed spectra where we indicate the spectral indices of both the inertial and kinetic ranges in Figure 8. This is essentially a limitation having multiple spacecraft at only a single length scale. Future multi-spacecraft plasma turbulence missions should be designed to sample multiple length scales simultaneously.
All spectra show a steepening of the spectra toward directions close to the mean magnetic field direction. Here we define the mean magnetic field to be the mean over the entire interval. A tendency toward a steepening in the parallel direction has been observed at larger scales [9] by varying investigating different orientations of the magnetic field and also with multi-spacecraft analysis at dissipation scales [42]. One of the predictions of a critically balanced cascade [80] is a spectral index relation for the magnetic field such that a Kolmogorov -5/3 spectrum in the perpendicular direction corresponds to a spectral index of -2 in the parallel direction in the inertial range. Extensions of this theory which assume a cascade of critically balanced kinetic Alfvén waves predict a parallel scaling of -5 and a perpendicular scaling of -7/3 [72,81,82]. Although the perpendicular scaling of the magnetic field is close to the -7/3 the parallel spectral index is far from -5, suggesting that this description is incomplete.
The presence of intermittency [83] or Landau damping [84] could modify the predictions of spectral indices to -8/3 in the perpendicular direction. Following the scaling relation k ∝ k 2/3 ⊥ this would correspond to a scaling of -7/2 in the parallel direction [83] which are more consistent with the observations reported here. Alternatively, the nature of the anisotropy in the spectral indices could be explained by a non-elliptic geometry of the power distribution without the need for critical balance [85].
An interesting feature shown here is that the density spectrum and the temperature spectra show similar evolution of the spectral index from ∼ −2.0 to ∼ −1.5. The similarity of the density and the thermal speeds is consistent with the interpretations of Gogoberidze et al. [73], although this comparison is not completely clear when investigating only the one dimensional PSD as a function of frequency shown in Figure 2. The extreme steepness in the velocity spectrum could be due to the ions becoming demagnetized and no longer participating in the cascade or it could be interpreted in terms of the Alfvén ratio which decreases rapidly and can be interpreted as a signature of a kinetic Alfvén wave [56,70,86].

Summary
In this study we have extended the MSR analysis technique to turbulent fields and scalars other than the magnetic field in the magnetosheath. The fields surveyed show both, power anisotropy, and anisotropy of spectral index. The anisotropy is the strongest in the velocity fluctuations which exhibits a steep power law at ion scales. The electron velocity (not shown) shows a much flatter spectra in the ion dissipation range of f −0.8 compared with the ion velocity spectra (f −3.08 ) indicating the ions are no longer magnetized whereas the electrons remain magnetized indicating differences in the ion and electron velocities and that the Hall effect is present. However, the electron dynamics is outside of the scope of this work but will be investigated in the future, with a smaller choice of inter-spacecraft distance. The compressible component of the turbulence appears to be less anisotropic for the magnetosheath. This is in contrast to some observations in the solar wind [22]. This suggests that the compressibility is different to that in the solar wind and likely to be due to compressible structures that are quasi aligned with the magnetic field direction.
The anisotropy of spectral index is investigated, finding that the spectra are steeper in the direction parallel to the mean magnetic field direction for all parameters measured. The vector fields show the largest variations, while the scalars show smaller variations. When compared to the 1D frequency spectra the spectral indices are similar except in the case of density which is most likely due to the small difference in the inertial range scaling and the transition range which comes between the inertial and dissipation ranges.
The thermal speed spectra and the density spectra show similar evolution of the spectral indices from parallel orientations to perpendicular ones suggesting that the prediction of Gogoberidze et al. [73] is correct. Finally the similarity of the density spectra with the thermal speed spectra may hint that particle heating and the energy cascade rate has a strong link to the compressibility of the plasma (e.g., [87,88]). There is also evidence that compressible slow waves and their interactions with Alfvén waves are simultaneously observed with heating signatures in proton velocity distribution functions [68]. These studies all highlight the importance of compressibility in understanding turbulent heating in a plasma.
There is also a hint of non-gyrotropy, although its origins and significance remain unclear. For this interval the anisotropy and non-gyrotropy appear to increase with one another, however more observations will be required to draw any conclusions from this case study. In the future the effect of the plasma beta, and the type of plasma should be investigated (i.e., fast/slow solar wind inner/outer magnetosheath etc.) which can all affect the level of anisotropy (e.g., [18,89]). These measurements are only possible due to the unique capabilities of MMS combining both multiple sampling points and high time resolution plasma measurements.
Here we investigate the ion dynamics, and will extend this work to electron dynamics when the stationarity of the fluctuations in the electron velocity is justified for the analysis. The structure of the fluctuations at ion scales are similar to one another irrespective of the parameters. Finally we have demonstrated that fully three dimensional observations of plasma parameters in the Earth's magnetosheath are possible with MMS, however one of the key limitations here is the ability to only study a single scale at once. Ideally all the relevant scales of the turbulence should be sampled simultaneously, this requirement should be a driver for the design of future spacecraft missions.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study can be found in the MMS science data archive https://lasp.colorado.edu/mms/sdc/public/.

AUTHOR CONTRIBUTIONS
OR anlyzed the data and drafted the manuscript. YN, ZV, and RN contributed to the interpretation of the analysis and general improvements in the manuscript. DG contributed to the interpretation of the FPI data and general improvements in the manuscript.

FUNDING
RN was supported by Austrian FWF projects I2016-N20. ZV was supported by the Austrian FWF projects P28764-N27.

APPENDIX Normalization
The result from the MSR technique is a power spectrum as a function of spacecraft frame frequency ω sc and wavevector k. For this example velocity units will be used: [(km/s) 2 (s rad)] (A1) To estimate the spectrum P MSR (k) an integration needs to be performed over ω sc , P MSR (k) = P MSR (ω sc , k)dω sc [(km/s) 2 ] (A2) The power spectrum can be expressed in cylindrical coordinates: To convert to a one dimensional spectrum the spectrum is averaged azimuthally: (A4) This is then folded along the line k = 0 (summed together and divided by 2) so there are only positive such that there are only positive components of k .
A dimension conversion to spectral density in the wavevector domain; where the prime denotes the change of units. A one dimensional cut of the spectrum is taken such that: Finally normalization to the fluctuation variance which assumes ergodicity of the signals reads: E MSR 1D (k ) = σ 2 P MSR Total (k ) P ′ MSR 1D (k ) [km 3 /(s 2 rad)] (A7) where σ 2 is the variance measured in the time domain: The spectrum in Equation (A7) normalized to the variance over the total power: