For our analysis, we use the dataset discussed by Alexandrova et al. (2012), when Cluster was in the free solar wind, i.e., not magnetically connected to the Earth’s bow shock. Details have been described also in Lacombe et al. (2017) and we recall here the main aspects. Magnetic field fluctuations are measured by the STAFF (Spatiotemporal Analysis of Field Fluctuation) instrument, composed of a waveform unit (SC) and a Spectral Analyzer (SA). Power spectra are computed on board in a magnetic field-aligned system of coordinates (MFA), based on the 4 s magnetic field measured by the FGM (Fluxgate Magnetometer) experiment. A selection of 112 spectra has been performed, retaining in each spectrum only measurements above three times the noise level in every direction x , y , and z (see Appendix in Lacombe et al., 2017). Each sample is a 10 min average of 150 individual 4 s spectral measurements. This provides spectra above 1 Hz up to typically 20–100 Hz, depending on the amplitude of the fluctuations in each interval. When converted into physical length scales, assuming the Taylor hypothesis ( k = 2 π f / V s w ), this leads to signals that cover the range between ∼ 2 d p and ∼ 0.5 d e (where d p and d e are the proton and electron inertial lengths, respectively), enabling then a good description of the sub‐ion regime from proton to electron scales.

The reference frame adopted (MFA) is such that B z is the component aligned with the mean magnetic field B 0 (relative to the 4 s interval during which an individual spectrum is calculated); B x is the component orthogonal to B z in the plane containing both the solar wind velocity V s w and the mean magnetic field B 0 ; and B y is the third orthogonal component. Note that a selection criterium is imposed on the angle θ B V , the angle between the local 4 s magnetic field, and the flow velocity; i.e., that θ B V is large enough to avoid a connection with the Earth bow shock during the sampled interval; θ B V in the dataset has an average value of ∼ 80 °. This implies that, for each spectrum, the mean magnetic field makes a big angle with respect to the sampling direction; moreover, we have checked that θ B V does not vary significantly during the 10 min over which spectra are averaged.

As a consequence, this procedure selects intervals in which Cluster observed highly oblique k-vectors and, to a good approximation, the component B x corresponds also to the sampling direction (radial) and is orthogonal to B 0 ; B y corresponds to the other perpendicular component; and B z is identified as the compressive component B ∥ . As already discussed in Lacombe et al. (2017), although the total trace power measured in situ is an invariant observable, the fact that the sampling occurs only in a preferred direction introduces a relative weight between B x and B y that is measurement dependent (Saur and Bieber, 1999). To take this into account, we have employed an analogous approach in the analysis of the simulations data, as described in the next section.

In situ observations are directly compared with numerical simulations performed with the hybrid-PIC code CAMELIA (Matthews, 1994; Franci et al., 2018a). Despite the fact that the hybrid model neglects the dynamics of electrons, it captures well the transition from fluid to kinetic regime around ion scales where electron effects do not play an important role. Hybrid simulations reproduce successfully many of the main properties of solar wind turbulence observed by spacecraft at sub‐ion scales (e.g., Perrone et al., 2013; Valentini et al., 2014; Franci et al., 2015a; Franci et al., 2015b; Franci et al., 2018b; Cerri et al., 2016; Cerri et al., 2017; Arzamasskiy et al., 2019). It is then a suitable tool to investigate the turbulent regime probed by STAFF/Cluster data. We use here 2D simulations—computationally more affordable than 3D—in order to explore the parameter space observed in situ; in particular, we focus on the effects associated with variations in the proton and electron plasma beta β p and β e . The restricted 2D geometry clearly cannot fully capture the richness of the turbulent phenomena (e.g., Howes, 2015) and in general, kinetic aspects related to the propagation of the fluctuations along the magnetic field are inhibited, like the presence of parallel propagating ion-scale waves and associated cyclotron resonances or the development of some kinetic instabilities (although some of their aspects can be still described also in 2D, e.g., Hellinger et al., 2015; Hellinger et al., 2017). On the other hand, in the case of the highly anisotropic solar wind turbulence, spectral properties can be captured efficiently (Franci et al., 2015a; Franci et al., 2015b). In particular, for the purpose of this work, Franci et al. (2016) have shown that 2D hybrid simulations are able to reproduce the ion-break scale behavior in different beta regimes observed in solar wind turbulence (Chen et al., 2014). Moreover, 3D hybrid simulations (Franci et al., 2018b) have confirmed the solidity of the reduced 2D results and the good agreement with in situ observations. We then exploit the good matching between simulations and in situ observations to characterize further the properties of kinetic plasma turbulence in the sub‐ion regime.

In order to make a direct comparison with sub‐ion spectra measured by Cluster, we have adopted a similar approach in the computation of spectra in the simulations. This means that numerical spectra are computed along the x direction only, to mimic the radial sampling occurring in the solar wind. This is obtained by integrating along y the Fourier spectrum P ( k x , y ) of each i magnetic field component: P i ( k x ) = ∫ ​ P i ( k x , y ) d y . (1) Therefore, also in the simulation, B x corresponds to the sampling direction, orthogonal to the out-of-plane magnetic field B z , and B y is the most energetic fluctuating component, being orthogonal to both B 0 and k = k x . With this approach and within the observational conditions previously described, we can perform a direct comparison of simulations and in situ data.

The numerical dataset used was originally presented in Franci et al. (2016) and is available online. It is constituted by a set of different 2048 2 2D simulations of decaying turbulence, corresponding to a physical simulation box size of 256 2 d i , except for the higher beta case, β p = 8 , where the size is 512 2 d i , and for different beta conditions covering the range of variations observed in situ, with β p = β e . Runs are initiated with random perpendicular Alfvénic fluctuations with vanishing cross-helicity and equipartition in magnetic and kinetic energies. The rms of the in-plane fluctuations is B r m s = 0.24 B 0 and the highest initially excited k-vector is k i n j d i = 0.2 ( B r m s = 0.48 B 0 and k i n j d i = 0.05 for the β p = 8 case). Spectra are computed at the maximum of the turbulent activity.

Bieber et al. (1996) and Saur and Bieber (1999) have investigated how different types of k-vectors distributions can generate a variable anisotropy in the observed magnetic field components, due to sampling effects. In the case of a gyrotropic 2D distribution of k-vectors, the ratio P y / P x is expected to coincide with the local slope γ of the spectrum P ( k ) ∼ k − γ . This applies well to solar wind observations in the physical range of interest here, as it can be appreciated in Figure 1, where the ratio P y / P x , shown in the bottom panels, is close to the spectral slope observed—typically in the range [−2.5,−3]—and appears roughly independent of the plasma beta. Interestingly, at smaller scales, when the magnetic spectrum steepens as approaching electron scales (Alexandrova et al., 2009), this is not associated with an increase in the perpendicular power ratio P y / P x (which on the contrary has a slight decrease); this does not correspond to the expectation for a quasi-2D spectrum according to the model and in fact, Lacombe et al. (2017) have interpreted this signature as the result of a more isotropic distribution of k-vectors close to electron scales.

To validate further this observational conclusion, we verify here the applicability of the Saur and Bieber model to sub‐ion scale turbulence. In the simulations, the spectrum is two-dimensional by construction and consistent with the axisymmetric initial conditions imposed in the x-y plane, it is also gyrotropic with respect to the out-of-plane magnetic field B z .

First, it is instructive to discuss spectra shown in Figure 3. These are power spectra of the perpendicular components B x (purple) and B y (orange) as a function of k x , assuming then a fixed direction of sampling. As expected, P y ( k x ) > P x ( k x ) ; on the other hand, their sum P ⊥ ( k x ) (solid black line) is statistically equivalent to the axisymmetric spectrum P ⊥ ( k ) = P y ( k ) + P x ( k ) . The difference is that when calculating the axisymmetric spectrum P ⊥ ( k ) , all perpendicular magnetic field directions have equal weight and one can assume that statistically P y ( k ) ∼ P x ( k ) ; as a consequence, the power associated with any individual perpendicular component corresponds to half of the total perpendicular power P ⊥ ( k ) / 2 ∼ P ⊥ ( k x ) / 2 (thin dashed black line). It is interesting to note that when sampling along a fixed direction (x), as it happens with spacecraft in the solar wind, none of the two measured spectra P y ( k x ) and P x ( k x ) is really representative of the power P ⊥ ( k ) / 2 of the gyrotropic description; instead, the component along the sampling ( B x ) is significantly reduced due to the solenoidal ∇ ⋅ B = 0 condition, while the orthogonal ( B y ) is amplified, in order to maintain the same total power P ⊥ ( k ) . This means that, in solar wind spectra like in Figure 1, neither P x nor P y is individually representative of the average power in a perpendicular B component: the individual measurements of P x or P y cannot be directly associated with it, but only their sum.

Bearing this in mind, Figure 4 shows the ratio of the power in the perpendicular components for the three simulations shown in Figure 2. The P y / P x ratio captures well the transition from MHD to a steeper spectrum at smaller scales; in all cases, the ratio, close to 5/3 at large scales, starts increasing in the vicinity of ion scales and reaches a maximum in the sub‐ion regime, where it is saturated close to ∼ 3 , in good agreement with the local spectral slope observed in the kinetic range, which is typically close to − 3 . At larger k, the ratio then decreases due to the noise. In the framework of the spectral anisotropy, Saur and Bieber model all this indicates a quasi-2D gyrotropic spectrum of the fluctuations, which corresponds well to the spectrum developed in these simulations. This confirms that the model is valid also at sub‐ion scales and reinforces the finding of Lacombe et al. (2017), where is found that solar wind spectra at kinetic scales are described well by a quasi-2D gyrotropic distribution.

There is another interesting indication suggested by Figure 4, namely, the fact that the P y / P x ratio in the sub‐ion range seems to depend on beta: consistent with this, the sub‐ion slope in Figure 2 is slightly steeper for small β p and shallower for larger β p . This behavior is already discussed in Franci et al. (2016) and is found in all simulations for the spectrum of the transverse fluctuations B ⊥ ; conversely, the spectrum of the parallel component B ∥ is almost independent of β p (see Figure 4 in Franci et al., 2016). We have then looked for a similar trend also in the in situ data. Figure 5 shows the histogram of the spectral slopes in the kinetic range for B ⊥ (top) and B ∥ (bottom), for larger (red) and smaller (black) total beta. Spectral slopes are calculated between 2 < k d p < 8 for β p < 1 and between 2 < k ρ p < 8 for β p > 1 , where a quite well-defined power-law scaling is observed. They are then separated into two groups defined by the total beta β < 2 and β > 2 . The mean of each histogram is indicated by the small vertical line ended with a diamond. For the parallel component (bottom panel), the distribution of the slopes is similar for both beta regimes and centered around a value of approximately − 2.65 ± 0.15 ; this is in good agreement with the simulations. For the dominant perpendicular component (top panel), we observe average values consistent with previous studies based on the total power δ B 2 = δ B ∥ 2 + δ B ⊥ 2 of the fluctuations (Alexandrova et al., 2009; Alexandrova et al., 2012; Chen et al., 2013a; Sahraoui et al., 2013). However, in the lower beta case (black), some slightly steeper slopes are observed for B ⊥ with respect to the high beta case, with an average of − 2.8 ± 0.15 with respect to − 2.7 ± 0.15 . We have checked that the difference in the histograms is statistically significant, thus suggesting some β-dependence in the spectral slope. A more detailed investigation is needed to fully identify the role of β p on the sub‐ion spectral slope and is beyond the scope of the present study. This behavior, however, agrees qualitatively with the simulations.

Moreover, a consequence of the behavior in Figure 5 is that while at high beta, δ B ∥ and δ B ⊥ have basically the same scaling, so that their ratio remains approximately constant in the sub‐ion range, at lower β, their slightly different scaling is expected to result in a slow increase of the δ B ∥ / δ B ⊥ ratio between ion and electron scales. These properties are related to the evolution of the magnetic compressibility of the fluctuations in the sub‐ion range, which is the main focus of the next section.

First, it is useful to go again from frequency to k-vector spectra: in Figure 8, frequencies are converted into k-vectors and normalized with respect to the proton inertial length d p .

We first identify two big categories such that both proton and electron betas are small, i.e., β p < 1 and β e < 1 , or both are large, i.e., β p > 1 and β e > 1 . We obtain an average total beta β ∼ 1 in the former and β ∼ 4 in the latter. The average spectrum of magnetic compressibility for each of the two families is shown in the top panel of Figure 8 as a function of k d p ; the thin dotted lines identify the standard deviation around the averages. In the high beta case (solid blue), the compressibility reaches a plateau after k d p = 1 and is saturated at an average level which is very close to isotropy (same power in P x , P y , and P z ), while in the low beta case (dashed red), C ∥ remains smaller and there is not a clear plateau at k d p > 1 .

The remaining spectra are further separated in two other families: the first with β p < 1 and β e > 1 and the second with β p > 1 and β e < 1 . In this case, the average total betas are very similar, β ∼ 1.9 ( β p ∼ 0.75 ) and β ∼ 2.0 ( β p ∼ 1.5 ), respectively, and fall in between the other two groups (small and large β). Consistent with this, the average spectrum of these two families, shown in orange and green in the bottom panel, has a level of compressibility at sub‐ion scales that is intermediate with respect to the other two curves. Moreover, they almost precisely fall on top of each other. All this suggests that not only is the total plasma beta a good parameter for ordering the level of compressibility generated at sub‐ion scales, but also this level is roughly independent of the individual weights of β p and β e , being their sum β = β p + β e the only relevant parameter.

This observational finding is in very good agreement with the expectation from the following relation: C ∥ = β p / 2 ( 1 + T e / T p ) 1 + β p ( 1 + T e / T p ) = β / 2 1 + β , (2) where T e and T p are the electron and proton temperatures.

Eq. 2 can be derived (Schekochihin et al., 2009; Boldyrev et al., 2013) under the assumption of low-frequency magnetic structures in pressure balance at scales where the ion velocity becomes negligible compared to the electron one, or equivalently, the Hall term J × B becomes dominant over the ideal MHD term − U × B . A special case is the regime of KAW, however, Eq. 2, which does not depend explicitly on k and thus on a specific dispersion relation, can be seen as a more general condition for highly oblique fluctuations in the sub‐ion range (e.g., ion-scale Alfvénic vortices, Jovanovic et al., 2020), under the assumptions described above (see e.g., Appendix C2 of Schekochihin et al., 2009).

To improve our analysis, we focus more in detail on the Cluster observations and compare them with numerical results. Note that, as in the simulations of Franci et al. (2016), it is only considered the case β p = β e ; we have made a selection of solar wind spectra with similar properties ( β p ∼ β e ∼ β / 2 ). These have then been divided in five subgroups as a function of β and averaged to obtain a mean C ∥ profile for each β-family. The selection results in 7, 13, 23, 9, and 1 spectra for β = 0.6 , 1,2,4,8 , respectively (only one spectrum fulfills the condition for high enough beta). Simulations with approximatively the same β p (and β) are considered for a direct comparison. In the following analysis, we want to identify the physical scale associated with the changes in the properties of the fluctuations and its possible connection to either the ion Larmor radius ρ p or the inertial length d p , as they are related by ρ p = β p d p .

The results of this comparison are shown in Figure 9, where scales are normalized to both d p (top) and ρ p (bottom). Left panels show spectra from in situ data and right panels result from simulations, where the colors encode the same range of β. Qualitatively, the global trend seen in the simulations matches well that of the observations. First, the level of magnetic compressibility reached sub‐ion scales increases monotonically with β, as expected. Second, we can identify a plateau phase beyond ion scales whose extension is gradually reduced as β decreases; for the smallest betas, the plateau disappears and is replaced by an almost monotonic increase of C ∥ all along the sub‐ion range, though with a shallower slope compared to that of the transition from the MHD range.

This seems to suggest a different behavior of the turbulent fluctuations populating the sub-ion cascade as a function of the beta. To investigate further this aspect, horizontal dotted lines in the right panels of Figure 9 show the theoretical prediction for the asymptotic level of C ∥ between ion and electron scales predicted by Eq. 2, with the same color scale. For simulations at large β, when a plateau is clearly observed, the level of magnetic compressibility also agrees well with the one predicted by the theory. In the low beta case, there is a larger discrepancy and the observed level of magnetic compressibility is larger than the constant level predicted by Eq. 2. The different behavior of the compressibility in low- and high beta regimes found in our simulations, together with the larger discrepancy with respect to the theoretical predictions observed at low beta, is also consistent with results from previous numerical studies (e.g., Cerri et al., 2016; Cerri et al., 2017; Grošelj et al., 2017).

The situation is somewhat different when comparing predictions to the in situ data; in this case, there is a slight difference between the KAW level and the observed one, and this is persistent at all β. In particular, at high beta, it is apparent that while Eq. 2 predicts compressibility that goes beyond 1/3 (for β → ∞ , we have C ∥ = 0.5 , so δ B ∥ = δ B ⊥ ), a condition well recovered in the simulations, in Cluster data C ∥ , does not go beyond component isotropy ( δ B ∥ = δ B ⊥ / 2 ; thus, C ∥ = 1 / 3 ). However, due to the low statistics in the data (just one spectrum has β ≳ 8 ), it is hard to draw a firm conclusion here.

Interestingly, from Figure 9, it seems that neither d p nor ρ p is able to fully capture and order the change in the spectrum of the magnetic compressibility for different betas; the saturation/plateau phase for low β spectra results more shifted toward high k-vectors compared to the high β ones when normalizing to d p , while the vice versa is observed when normalizing to ρ p . This suggests that the behavior can be better captured by an intermediate scale between the two. For this reason, in Figure 10, we have normalized spectra on a mixed scale d p ρ p . Note that such a scale, proportional to d p β p 1 / 4 , was found to describe well the behavior of the ion-break scale in magnetic field spectra in the range β p ∼ 1 by Franci et al. (2016), and, although not shown, to describe the variation of the break of the parallel magnetic field spectrum at all betas; this then motivated our choice. When such a mixed scale is used (top right panel), all cases follow the same trend: they grow until they reach k d p ρ p ∼ 2 and then start flattening, the saturation level depending on the beta. In situ observations (top left panel) seem to follow the same trend, confirming that such an intermediate scale is a good candidate for controlling the variation of the magnetic compressibility spectrum at ion scales.

It is then reasonable to use such a k-vector normalization to better evaluate the agreement with Eq. 2. In the bottom panels of the same figure C ∥ * , spectra are then normalized to the theoretical prediction for C ∥ . In simulations, as already pointed out, cases with β > 1 display a good agreement with the sub‐ion compressibility level predicted by the theory; as a consequence, when normalized to d p ρ p , all spectra collapse on top of each other all along ion and sub‐ion scales. A worse agreement is observed at β ≤ 1 when simulations display a slightly higher compressibility level than predicted. Quite differently, the ratio between the in situ observations and the theoretical C ∥ is always below one and around 0.7–0.8 for all β groups in the sub‐ion range (see also Figure 10 of Lacombe et al., 2017). This behavior is consistent with the results of Pitňa et al. (2019) based on observations from the wind spacecraft, who find on average C ∥ ∼ 0.9 , without making a distinction among beta regimes and with most of the data displaying a slightly smaller magnetic compressibility than the prediction. Our study confirms this scenario and suggests that the same trend is followed for all spectra, almost independently of the plasma beta. A ratio smaller than one and close to ∼ 0.75 is also consistent with similar observational results of the plasma compressibility and based on the ratio between density and perpendicular magnetic fluctuations predicted by linear theory (Chen et al., 2013a; Pitňa et al., 2019). This was interpreted by Chen et al. (2013a) as a consequence of the nonlinear behavior of the solar wind fluctuations in the sub‐ion range, in agreement with simulations of strong KAW-turbulence (Boldyrev et al., 2013). On the other hand, for the magnetic compressibility, our fully nonlinear simulations of sub‐ion turbulence do not recover the same effect seen in situ, as C ∥ * ≳ 1 . Other reasons could explain such a discrepancy, e.g., the effect of some electron Landau damping on the fluctuations observed in situ (Howes et al., 2011; Passot and Sulem, 2015; Schreiner and Saur, 2017) and not captured by the hybrid model. In order to answer these questions, a more detailed study of the polarization properties of the fluctuations in our simulations is in preparation.

Finally, note that the increase in C ∥ observed at higher k in the in situ data could be related to a further change in the properties of the fluctuations as they approach electron scales; as discussed in Lacombe et al. (2017), this also coincides with a change in the estimated spectral anisotropy. For example, Chen and Boldyrev (2017) have suggested that the increase in the magnetic compressibility beyond the sub‐ion range could be related to electron inertia corrections to Eq. 2. This effect is then not captured by the hybrid model and we cannot compare any more the observations with the simulations in this range. It is however interesting to note that while the further increase of compressibility at electron scales is predicted for β e ≲ 1 (Chen and Boldyrev, 2017; Passot et al., 2017), in the intervals measured by Cluster, it seems to be observed for all beta ranges for k d p ≳ 10 ( k d e ≳ 1 / 4 ). Moreover, it is also interesting to note that spectra for all betas reach isotropy at roughly k ρ p ∼ 20 , corresponding on average to k ρ e ∼ 0.5 .