Magnetic field turbulence in the solar wind at sub-ion scales: in situ observations and numerical simulations

We investigate the transition of the solar wind turbulent cascade from MHD to sub-ion range by means of a detail comparison between in situ observations and hybrid numerical simulations. In particular we focus on the properties of the magnetic field and its component anisotropy in Cluster measurements and hybrid 2D simulations. First, we address the angular distribution of wave-vectors in the kinetic range between ion and electron scales by studying the variance anisotropy of the magnetic field components. When taking into account the single-direction sampling performed by spacecraft in the solar wind, the main properties of the fluctuations observed in situ are also recovered in our numerical description. This result confirms that solar wind turbulence in the sub-ion range is characterized by a quasi-2D gyrotropic distribution of k-vectors around the mean field. We then consider the magnetic compressibility associated with the turbulent cascade and its evolution from large-MHD to sub-ion scales. The ratio of field-aligned to perpendicular fluctuations, typically low in the MHD inertial range, increases significantly when crossing ion scales and its value in the sub-ion range is a function of the total plasma beta only, as expected from theoretical predictions, with higher magnetic compressibility for higher beta. Moreover, we observe that this increase has a gradual trend from low to high beta values in the in situ data; this behaviour is well captured by the numerical simulations. The level of magnetic field compressibility that is observed in situ and in the simulations is in fairly good agreement with theoretical predictions, especially at high beta, suggesting that in the kinetic range explored the turbulence is supported by low-frequency and highly-oblique fluctuations in pressure balance, like kinetic Alfv\'en waves or other slowly evolving coherent structures.


INTRODUCTION
The solar wind constitutes a unique laboratory for plasma turbulence (Bruno and Carbone, 2013). In the last decade increasing interest has raised towards the small-scale behaviour of the turbulent cascade, i.e. beyond the breakdown of the fluid/MHD description that takes place at ion scales. Spacecraft observations of solar wind and near-Earth plasmas provide unique measurements of the turbulent fluctuations at scales comparable and smaller than the typical particle scales, the Larmor radius ρ (see Appendix for definition of physical quantities used) and the inertial length d (e.g. Alexandrova et al., 2009;Sahraoui et al., 2010;Alexandrova et al., 2012;Chen et al., 2013a). However, the physical processes governing the energy cascade at kinetic scales and those responsible for its final dissipation are not well understood yet.
What is well established is that in the transition from MHD to kinetic regime, plasma turbulence modifies its characteristics. Observational and numerical studies over the last few years have highlighted the main differences between large and small scale properties of solar wind fluctuations (e.g., Chen, 2016;Cerri et al., 2019). The magnetic field spectrum typically steepens when approaching ion scales, leading at sub-ion scales (between ion and electron typical scales) to a power law with spectral index close to −2.8 (Alexandrova et al., 2009(Alexandrova et al., , 2012Kiyani et al., 2009;Chen et al., 2010;Sahraoui et al., 2013), steeper than Kolmogorov -5/3, but also than the theoretical prediction −7/3 from EMHD (Biskamp et al., 1996) and KAW/whistler turbulence (Schekochihin et al., 2009;Boldyrev et al., 2013). The origin of such a spectral slope is still unknown and it has been proposed that it could be related to intermittency corrections (Boldyrev and Perez, 2012;Landi et al., 2019), magnetic reconnection (Loureiro and Boldyrev, 2017;Mallet et al., 2017;Cerri et al., 2018), Landau damping (Howes et al., 2008;Schreiner and Saur, 2017) and/or the role of the non-linearity parameter (Passot and Sulem, 2015;Sulem et al., 2016).
The change in the magnetic field spectrum is accompanied by a rapid decrease in the power of ion velocity fluctuations (Šafránková et al., 2013;Stawarz et al., 2016) and the onset of the non-ideal terms in the Ohm's law which governs the electric field associated to the turbulent fluctuations; as a consequence the electric field spectrum becomes shallower at sub-ion scales (Franci et al., 2015a;Matteini et al., 2017). In this framework the electric current (mostly carried by electrons) plays a major role, coupling directly with the magnetic field in the cascade and likely affecting the energy cascade rate via the Hall term (Hellinger et al., 2018;Papini et al., 2019). All these properties depend further on the plasma beta (β = 8πnk B T /B 2 ), which controls, among other things, the scale at which the magnetic field spectrum breaks (Chen et al., 2014;Franci et al., 2016).
One of the most significant differences with respect to the turbulent regime observed at large scales however is the role of compressive effects. While in the inertial range fluctuations show a low level of both plasma and magnetic field compressibility, hence can be reasonably well described by incompressible MHD, at sub-ion scales density and magnetic field intensity fluctuations become significant and comparable to transverse ones (Alexandrova et al., 2008;Sahraoui et al., 2010;Salem et al., 2012;Kiyani et al., 2013;Chen et al., 2012b;Perrone et al., 2017), in agreement with simulations (Franci et al., 2015b;Cerri et al., 2017). It is believed that this is related to a change in the properties of the turbulent fluctuations, which become intrinsically compressive at small scales. It is then by studying in detail their properties that it is possible to shed light on the nature of the fluctuations which support the cascade at kinetic scales (Chen et al., 2013b;Pitňa et al., 2019;Grošelj et al., 2019;Alexandrova et al., 2020).
Another important aspect of solar wind turbulence is its spectral anisotropy (Horbury et al., 2008;Wicks et al., 2010;Chen et al., 2010;Roberts et al., 2017). Studies about the shape of turbulent eddies, both at MHD (Chen et al., 2012a;Verdini et al., 2018Verdini et al., , 2019 and kinetic scales (Wang et al., 2020), reveal the presence of a 3D anisotropy in the structures when described in terms of a local frame. On the other hand, when the analysis is made in a global frame (without tracking the local orientation of the structures), the 3D anisotropy is not captured, and the k-vectors of the fluctuations show a statistical quasi-2D distribution around the magnetic field (Matthaeus et al., 1990;Dasso et al., 2005;Osman and Horbury, 2006). In this work we address this latter aspect and we investigate the distribution of the k-vectors with respect to the ambient magnetic field at kinetic scales by using the magnetic field variance anisotropy (i.e. the ratio of magnetic field fluctuations in different components). Saur and Bieber (1999) have shown that, also in single spacecraft observations, is possible to characterize the 3D k-vector distribution by using variance anisotropy. When the sampling occurs only along a preferential direction, like in typical solar wind observations, their model predicts various possible kinds of variance anisotropy as a function of the underlying k-spectrum. In particular, assuming a quasi-2D gyrotropic distribution of k-vectors (axisymmetric with respect to the magnetic field), the ratio of the power in the two perpendicular magnetic field components is directly related to the local slope of the spectrum -which is assumed to have the same form for all components and a slope independent of the scale within a given regime. Since both quantities, spectral slope and perpendicular power ratio, can be easily measured in situ, the Saur & Bieber model constitutes an useful and simple tool to investigate underlying spectral anisotropies. Despite the model was originally developed for MHD scale fluctuations, it basically corresponds to a geometrical description built on the divergence-less condition for B, so it can be applied to any kind of regimes, including the low-frequency turbulence expected at sub-ion scales (Turner et al., 2011). In the work of Lacombe et al. (2017), we investigated the k-vector distribution at sub-ion scales using the technique by Saur and Bieber (1999). Based on the comparison with the predictions, we concluded that the distribution of the k-vectors in the sub-ion range of solar wind turbulence is consistent with a quasi-2D gyrotropic spectrum, then approaching a more isotropic shape when reaching electron scales . However, such an application has not been benchmarked by kinetic numerical studies yet.
The aim of this work is then to focus on the spectral anisotropy properties and magnetic compressibility at small scales, by exploiting the detailed comparison of in situ observations and highresolution kinetic numerical simulations. The paper is organised as follow: in Sec. 2 we introduce the spacecraft and numerical dataset used and in Sec. 3 we describe their spectral properties. In Sec. 4 we discuss the spectral anisotropy at sub-ion scales and test, for the first time, the Saur and Bieber model in numerical kinetic simulations; in Sec. 5 we address properties of the magnetic compressibility and its dependence on the plasma beta. Finally, in Sec. 6 we discuss our conclusions and the implications of our finding for the interpretation of solar wind observations and simulations.

DATA AND SIMULATIONS
In this study we compare properties of magnetic fluctuations measured in situ by the Cluster spacecraft with numerical results obtained by means of 2D hybrid particle-in-cell (PIC) simulations.

Cluster STAFF spectra
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 (Spatio-Temporal Analysis of Field Fluctuation) instrument, composed by a wave form unit (SC) and a Spectral Analizer (SA). Power spectra are computed on board in a magnetic field aligned system of coordinates (MFA), based on the 4s magnetic field measured by the FGM (Fluxgate Magnetometer) experiment. A selection of 112 spectra has been performed, retaining in each spectrum only measurements above 3 times Figure 1. Cluster STAFF spectra for different intervals with β = 0.8, 1.5, 4 from left to right. Colors encode the magnetic field components B x (black dashed), B y (blue), and B z (red). The region highlighted in yellow corresponds to the the sub-ion range investigated in this study. Bottom panels show the ratio of the power in the perpendicular components P y (k)/P x (k) (black diamonds) and the value γ of the local slope of the total spectrum P (k) ∼ k −γ (red stars). The average values of P y (k)/P x (k) and γ in the sub-ion range are also shown as horizontal dashed lines. Figure 2. Magnetic field spectra from hybrid simulations for different beta regimes (β p = 0.125, 1, 10 and β e = β p ). Components are encoded as in Figure 1 and the coloured region indicates the sub-ion range that can be directly compared with the analogous region in the observations. the noise level in every direction x, y and z (see Appendix in Lacombe et al., 2017). Each sample is a 10 minute average of 150 individual 4s spectral measurements. This provides spectra above 1Hz up to typically 20-100Hz, depending on the amplitude of the fluctuations in each interval. When converted into physical length scales, assuming the Taylor hypothesis (k = 2πf /V sw ), this leads to signals that cover the range between ∼ 2d p and ∼ 0.5d 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 4s 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 sw 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 θ BV , the angle between the local 4s magnetic field and the flow velocity, i.e., that θ BV is large enough to avoid a connection with the Earth bow shock during the sampled interval; θ BV in the dataset has an average value of ∼ 80 degrees. 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 θ BV does not vary significantly during the 10 minutes 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.

Hybrid 2-D numerical simulations
In situ observations are directly compared with numerical simulations performed with the hybrid-PIC code CAMELIA (Matthews, 1994;Franci et al., 2018a). The hybrid model captures well the transition from fluid to kinetic regime around ion scales. Moreover, it reproduces successfully many of the main properties of solar wind turbulence observed by spacecraft at sub-ion scales (Franci et al., 2015b,a). 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 to variations in the proton and electron plasma beta β p and β e . Franci et al. (2016) have shown that 2D hybrid simulations are able to capture the ion-break scale behaviour in different beta regimes observed in solar wind turbulence (Chen et al., 2014). We then exploit the good matching between the simulations and in situ observations to characterise further the properties of kinetic plasma turbulence in the sub-ion regime. On the other hand, 3D hybrid simulations (Franci et al., 2018b) have confirmed the solidity of the reduced 2D results and the good agreement with in situ observations.
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: Therefore, also in the simulation B x corresponds to the sampling direction, orthogonal to the outof-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 2-D simulations of decaying turbulence with different beta conditions covering the range of variations observed in situ, and β p = β e ; runs are initiated with random perpendicular Alfvénic fluctuations with vanishing cross-helicity and equipartition in magnetic and kinetic energies. Spectra are computed at the maximum of the turbulent activity. Observations cover ion and electron scales, with a transition accompanied by a slope change around kd p ∼ 10 . In this work, we focus on the sub-ion regime highlighted in yellow in the panels, where electron physics effects can be neglected (at least for spectral properties) and a well-defined slope close to −2.8 can be observed (Alexandrova et al., 2012). The three cases, corresponding to different total beta β regimes [0.5, 1.5, 4], show a similar qualitative behaviour: as expected, the spectrum P y of the perpendicular B y component (blue) is always the most energetic. The power in the other perpendicular component P x (black dashed) is always slightly smaller, however, its ratio with P y is roughly independent of beta and close to the local spectral slope (bottom panels); this is related to the 3-D distribution of k-vectors  and will be discussed more in detail in Sec.4.

IN SITU DATA ANALYSIS AND SIMULATION RESULTS
On the other hand, the power P z of the field aligned component B z (red) is typically less energetic than P y , however, its weight is highly variable with beta: P z is smaller than P x for β < 1, comparable to P x for β ∼ 1, and larger the P x for β > 1. This obviously results in a variable magnetic compressibility associated to the fluctuations and its functional dependence on beta is the subject of Sec.5 Figure 2 shows an analogous selection from numerical simulations; note that in the simulations β e = β p . In this case the regime reproduced in the simulation box includes the MHD inertial range and its transition to a sub-ion cascade at smaller scales. The yellow area highlights the region of the Figure 3. Reduced spectra of the fluctuations of the magnetic field components B y and B x defined with respect to a fixed sampling direction k x for a simulation with β p = 0.5. The thick solid black line corresponds to the total perpendicular power P ⊥ (k x ); the dashed line shows P ⊥ (k x )/2, also corresponding to the average power in any perpendicular magnetic field component in the axisymmetric case. spectra − roughly a decade between kd p ∼ 1 and kd p ∼ 10 − that can be directly compared with the in situ data. In this region, the qualitative behaviour of the spectra is similar to Figure 1: B y (blu) is always dominant, B x (black) contributes for a constant fraction of it and is roughly the same at all betas, while B z (red) varies significantly in the panels and becomes comparable to B y for large betas. This confirms that our method of computing spectra in the simulations mimicking satellite observations really captures the main aspects of in situ measurements and can then be exploited to investigate further the properties of the turbulent cascade. 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 Fig. 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 to an increase in the perpendicular power ratio P y /P x (which on the contrary has a slight decrease); this doesn't correspond to the expectation for a quasi-2D spectrum according to the model and in fact Lacombe et al. (2017) has interpreted this signature as the result of a more isotropic distribution of k-vectors close to electron scales.

Perpendicular components ratio
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 to 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 . Hybrid simulations: spectra of the ratio of the perpendicular magnetic components P y (k x )/P x (k x ) and corresponding to the local spectral slope. Different colors encode different β p : 0.125 (cyan), 1 (black) and 8 (red). The horizontal dashed lines show reference spectral slopes observed in the simulations at kd p < 1 (-5/3) and at sub-ion scales kd p > 1.
(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 Fig. 1, neither P x nor P y are individually representative of the average power in a perpendicular B component: the individual measurements of P x or P y cannot be directly associated to 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 well captures 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 reach a maximum in the sub-ion regime, where the ratio saturates 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 Saur and Bieber model for spectral anisotropy 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.

Beta dependence
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 behaviour 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 < kd 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 in 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 centred 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., , 2012Sahraoui et al., 2013;Chen et al., 2013a). 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. This behaviour agrees qualitatively with the simulations; however 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. What can be however pointed out is that a consequence of the behaviour in Figure 5 is that while at high beta, δB and δB ⊥ are observed to have almost 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 main focus of the next section.

MAGNETIC COMPRESSIBILITY
We now investigate the role of the third magnetic field component B z , which is aligned with the local (at 4s) magnetic field B 0 . In particular, we focus on the the magnetic compressibility C = δB 2 /δB 2 , where δB 2 = δB 2 +δB 2 ⊥ and its implication for the nature of the cascade at these scales. Note that in this case, the measurement of B z is not affected by Figure 6. Examples of Cluster FGM (black) and STAFF (red) spectra of magnetic compressibility C as a function of the frequency measured in spacecraft frame.
the sampling direction (provided that this is orthogonal to B 0 to a good approximation) and since we use the total perpendicular power P ⊥ , the caution discussed in Sec. 4 is not needed here. Figure 6 shows C for three intervals of different total β=1,3,4 (β p = 0.3, 1.4, 2.5) as measured from STAFF (red). For these three cases we also show the spectrum of the magnetic field compressibility as measured at lower frequencies (corresponding to physical scales larger than d p ) by the fluxgate magnetometer onboard Cluster (FGM, black). Note that FGM spectra are linearly interpolated between 0.14 and 0.4Hz to remove artefacts due to spacecraft spin (0.25 Hz). There is a good matching between the two independent measurements at f ∼ 1Hz and where data points from both instruments are available for a more extended range there is also a quite satisfactory overlap between them. The overall behaviour agrees well with the expected picture: at large scale, in the MHD inertial range, the level of compressibility is lower, typically C 0.1 (e.g., Horbury and Balogh, 2001;Smith et al., 2006) and starts to increase as approaching ion scales (Hamilton et al., 2008;Alexandrova et al., 2008;Salem et al., 2012;Kiyani et al., 2013), reaching sometimes variance isotropy (indicated by the dashed horizontal line) in the sub-ion range, where the compressibility seems to saturate. As already shown by Figure 7. STAFF spectra of magnetic compressibility C as a function of the frequency measured in spacecraft frame. Different colours and lines identify different groups of intervals with given β. Lacombe et al. (2017), the level of magnetic compressibility developed at small scales is larger for high beta than for small beta. Since we focus on the behaviour at sub-ion scales, in the following we restrict our analysis to STAFF measurements only.
To highlight further the β-dependence of the magnetic compressibility, Figure 7 shows C for a selection of spectra with different β, increasing from red to purple. There is a continuous transition from lower to higher magnetic compressibility as a function of beta, in agreement with linear theory expectations (e.g. Podesta and TenBarge, 2012). Moreover, at high beta it seems that the fluctuations reach an asymptotic δB 2 /δB 2 ratio, leading to an extended plateau in the spectrum, while at the lowest beta a plateau cannot be clearly identified. We now want to identify more in detail what process and length scale control the level of C and in solar wind data.

Beta dependence and theoretical predictions
First it is useful to go again from frequency to kvector spectra: in Figure 8 frequencies are converted into k-vectors and normalised 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 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 kd 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 kd p = 1 and saturates 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 kd p > 1.
The remaining spectra are further separated in two other families: the first with β p < 1 and β e > 1, the second β 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 the total plasma beta is a good parameter for ordering the level of compressibility generated at sub-ion scales, but also this level is roughly independent on 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 relation below: 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 lowfrequency magnetic structures in pressure balance at scales where the ion velocity becomes negligible compared to the electron one, or equivalently Figure 8. Cluster average spectra of magnetic compressibility for intervals with β e and β p < 1 (red) and β e and β p > 1 (blue); in the bottom panel, intermediate values with similar average β but with β p < 1, β e > 1, and β p > 1, β e < 1 are also shown in green and orange respectively. Thin dotted lines in the upper panel show the one-sigma dispersion of the data.
the Hall term J × B becomes dominant over the ideal MHD term −U × B. A special case is the regime of kinetic Alfvén waves (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).

Comparison with simulations
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 5 sub-groups as a function of β and averaged to obtain a mean C Figure 9. Top panels: (Left) Cluster spectra of magnetic compressibility for intervals binned on different β, encoded in different styles and colours. Only cases with β p ∼ β e have been retained. The horizontal dashed line denotes energy equipartition between components (i.e. isotropy). (Right) Spectra of magnetic compressibility for simulations with different β p = β e , shown with same style as left panel. The increase of C for kd i 8 is due to numerical noise. Bottom: Same as top panels, but with k-vectors normalised with respect to the ion gyro-radius ρ p . Horizontal dotted lines, coloured according to their β, are the theoretical prediction of C from Eq.(2).
profile for each β-family. The selection results in 7, 13, 23, 9 and 1 spectra for β = 0.6, 1, 2, 4, 8 respectively (only 1 spectrum fulfils 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 to 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: The results of this comparison are shown in Figure 9, where scales are normalised to both d p (top) and ρ p (bottom). Left panels show spectra from in situ data and right panels results 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 at 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 behaviour 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 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 behaviour 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, are also consistent with results from previous numerical studies (e.g. Cerri et al., 2016Cerri et al., , 2017Groš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 a 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 1 spectrum has β 8) it's hard to draw a firm conclusion here.
Interestingly, from Figure 9 it seems that nor d p nor ρ p are 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 towards high k-vectors compared to the high β ones when normalising to d p , while the vice-versa is observed when normalising to ρ p . This suggests that the behaviour can be better captured by an intermediate scale between the two. For this reason, in Figure 10 we have normalised spectra to a mixed scale d p ρ p .
Note that such a scale, proportional to d p β 1/4 p , was found to describe well the behaviour of ion-break scale in magnetic field spectra in the range β p ∼ 1 by Franci et al. (2016), and, although not shown, to describes 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 normalised 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 normalised 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 1 and around 0.7-0.8 for all β groups in the sub-ion range (see also Figure 10 of Lacombe et al., 2017). This behaviour 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 Figure 10. Same spectra as in Figure 9, but with k-vectors normalised to the mixed scale d p ρ p ; in the bottom panels C is normalised to the theoretical prediction by Eq.2. same trend is followed for all spectra, almost independently of the plasma beta. A ratio smaller than 1 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 non-linear behaviour 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 non-linear 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 polarisation 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 kd p 10 (kd e 1/4). Moreover, it's also interesting to note that spectra for all betas reach isotropy at roughly kρ p ∼ 20, corresponding on average to kρ e ∼ 0.5.

CONCLUSION
In summary, we have discussed properties of magnetic field spectra of turbulent fluctuations in the sub-ion regime and their main dependence on the plasma beta. We have carried out a detailed comparison between in situ Cluster magnetic field observations in the frequency range f (Hz)= [1,200], corresponding to scales typically between d p < l < d e , and high-resolution 2D hybrid simulations.
First we investigated the spectral anisotropy of magnetic fluctuations at sub-ion scales. Our simulations confirm that the model of Saur and Bieber (1999), originally developed for MHD range fluctuations, is valid also at kinetic scales; by applying the model to the numerical spectra obtained mimicking the sampling along a fixed direction made by spacecraft, we were able to successfully capture original spectral properties as well as their variation with β. This then reinforces the finding of Lacombe et al. (2017) who applied the Saur and Bieber model to kinetic-scale observations for the first time and concluded that fluctuations of the solar wind spectrum in the sub-ion range are quasi-2D and gyrotropic. Moreover, we have shown that the component anisotropy measured in situ -leading to an apparent non-gyrotropic spectrum from an original gyrotropic one (see also Turner et al., 2011) is a direct consequence of the solenoidal condition of the magnetic field and the sampling procedure. This is not an effect related to the Doppler-shift of k-vectors swept through the spacecraft by the fast plasma flow and in fact, we were able to reproduce it in simulations just imposing a fixed sampling direction.
Note that our result about the global 2D-symmetry of the k-vectors around the magnetic field is not inconsistent with studies addressing the local shape of the eddies and suggesting the presence of a 3D anisotropy (e.g. Chen et al., 2012a;Verdini and Grappin, 2015;Verdini et al., 2018Verdini et al., , 2019Wang et al., 2020). In our approach we do not consider the specific orientation of the turbulent structures in the plane perpendicular to B, and it is reasonable to expect that their local 3D anisotropy is then lost. In other words, despite the possible presence of a 3D anisotropy of the turbulent eddies their k-vectors can be oriented isotropically around B, leading then -in a frame like the one used here -to the 2D spectrum found in the Cluster observations. This does not exclude that some aspects of the 3D anisotropy could be still captured using also a global approach, however, our study suggests that in this case one has to take carefully into account the effects of the apparent component anisotropy introduced by the sampling (Saur and Bieber, 1999, see also Figure 3 in this work).
For the magnetic compressibility C , we have confirmed that it has a strong dependence on the plasma beta (e.g., Alexandrova et al., 2008;TenBarge et al., 2012;Lacombe et al., 2017). In particular we have shown that in Cluster observations C depends on the total beta β only (Fig. 8), as expected for lowfrequency pressure-balanced fluctuations at highly oblique propagation (e.g., KAW). In the β range explored we find a good qualitative agreement between the trend observed in the data and in the simulations. The compressibility is observed to increase as a function of β, leading to a plateau at sub-ion scales for high betas and in good agreement with the prediction by Eq. (2). At low beta, a full developed plateau is not observed beyond ion scales and the compressibility continue to slowly increase along sub-ion scales, both in observations and simulations (see also Grošelj et al., 2019). There is, however, a difference in the asymptotic level of compressibility reached at high β in data and our simulations; in the former, fluctuations seem not to exceed component isotropy (C = 1/3), while in the latter they approach C = 0.5, which is the limiting value predicted by Eq. (2). This aspect deserves to be explored in future studies, extending the range of β explored, to then establish if the asymptotic condition observed in simulations and predicted by the theory, which implies same power in the parallel component as in the sum of the perpendicular ones, can be also observed in situ for high enough β intervals. As a consequence of the behaviour just described, there is a different quantitative agreement of the magnetic compressibility observed in situ and in simulations, with the theoretical prediction by Eq. (2). In simulations there is very good matching with the predicted level at higher beta, but an excess of C at low beta; this effect was already observed in Cerri et al. (2017) and is confirmed here on a large range of β. On the other hand, in solar wind observations, the ratio is always lower than 1 (smaller compressibility than predicted by the theory), and close to ∼ 0.75 for all β, in agreement with similar studies on the plasma compressibility (Chen et al., 2013a;Pitňa et al., 2019).
Our analysis also suggests that the increase in the compressibility at ion scales is controlled by an intermediate scale between the Larmor radius ρ p and the proton inertial length d p (Fig. 9). For simulations this was already anticipated in Franci et al. (2016), and we could identify it as related to d p ρ p , thus proportional to d p β 1/4 p (Fig. 10). Such a scaling with β p also corresponds to the scaling observed for the spectral ion-break in the range β p ∼ 1. However, it is worth to highlight that both observations (Chen et al., 2014) and our simulations (Franci et al., 2016) show that the spectral ion-break scale follows the largest of ρ p and d p depending on the beta, so that the correction term proportional to d p β 1/4 p identified in Franci et al. (2016) is important only around β p ∼ 1. On the other hand, the present study indicates that a scale proportional to d p ρ p orders well the spectra of compressibility at all betas, for both in situ data and simulations, suggesting that such a mixed scale controls the transition in the nature of the fluctuations from MHD to sub-ion range (see also the monotonic scaling with β p of the ion-break in the parallel magnetic field spectrum shown in Fig.4 of Franci et al., 2016). This may imply that the two changes of regimethe steepening of the magnetic spectrum and the increase in the compressibility -can occur at different scales for more extreme β values. In particular, we expect the spectral break to occur at a larger scale with respect to the plateau in the compressibility when β p 1 or β p 1, as in these cases d p ρ p is always smaller than the largest between ρ p and d p . A more detailed analysis on this aspect will be the subject of a future study, as well as the possible implications of this behaviour for fluctuations in the inner Heliosphere, where the plasma beta is typically lower than at 1AU, which can be observed by the Parker Solar Probe and Solar Obiter.

APPENDIX: SYMBOL DEFINITIONS AND NORMALIZED UNITS
The subscripts ⊥ and refer to the direction with respect to the ambient magnetic field B 0 and p and e denote respectively protons and electrons. All equations are expressed in the c.g.s. unity system. n and T denote the number density and the temperature of a species (we assume also n p = n e = n). β e,p = 8πnk B T e,p /B 2 0 are the electron and proton betas, and β = β p + β e is the total plasma beta; here k B is the Boltzmann constant. For each species of mass m and charge q, the inertial length d is defined a c/ω p , where ω p = (4πnq 2 /m) 1/2 is the plasma frequency, and the Larmor radius ρ is defined as v th /Ω c where v th it the thermal speed of each species and Ω c = q p B o /mc is the cyclotron frequency. V sw is the solar wind speed and f the frequency of the fluctuations measured by the spacecraft; k denotes the module of the wave vector k. presentation in the manuscript. All authors revised the manuscript before submission.

FUNDING
This work was supported by the Programme National PNST of CNRS/INSU co-funded by CNES. It has also been funded by Fondazione Cassa di Risparmio di Firenze through the project HY-PERCRHEL. We acknowledge PRACE for awarding us access to resource Cartesius based in the Netherlands at SURFsara through the DECI-13 (Distributed European Computing Initiative) call (project "HybTurb3D"), and CINECA for the availability of high performance computing resources and support under the ISCRA initiative (grants HP10C877C4 and HP10BUUOJM) and the program Accordo Quadro INAF-CINECA 2017-2019 (grants C4A26 and C3A22a). L.F. was supported by Fondazione Cassa di Risparmio di Firenze, through the project "Giovani Ricercatori Protagonisti", and by the UK Science and Technology Facilities Council (STFC) grants ST/P000622/1 and ST/T00018X/1. PH acknowledges grant 18-08861S of the Czech Science Foundation. OA and CL are supported by the French Centre National d'Etude Spatiales (CNES).