Characterizing current structures in 3D hybrid-kinetic simulations of plasma turbulence

In space and astrophysical plasmas turbulence leads to the development of coherent structures characterized by a strong current density and important magnetic shears. Using hybrid-kinetic simulations of turbulence (3D with different energy injection scales) we investigate the development of these coherent structures and characterize their shape. First, we present different methods to estimate the overall shape of the 3D structure using local measurements, foreseeing application on satellite data. Then we study the local magnetic configuration inside and outside current peak regions, comparing the statistics in the two cases. Last, we compare the statistical properties of the local configuration obtained in simulations with the ones obtained analyzing an MMS dataset having similar plasma parameters. Thanks to our analysis, 1) we validate the possibility to study the overall shape of 3D structures using local methods, 2) we provide an overview of local magnetic configuration emerging in different turbulent regimes, 3) we show that our 3D-3V simulations can reproduce the structures that emerge in MMS data for the periods studied by Phan et al. (2018), Stawarz et al. (2019).


Introduction
Turbulent, magnetised plasmas permeate a wide range of space and astrophysical environments, and plasma turbulence naturally develops coherent structures characterized by high current density and strong magnetic shear. These features are indeed present in practically any turbulence simulation employing the most disparate plasma models and regimes (e.g., Zhdankin et al. 2013;Passot et al. 2014;Navarro et al. 2016;Zhdankin et al. 2017;Cerri et al. 2019;Comisso and Sironi 2019, and references therein), as well as routinely observed via in-situ measurements in space plasmas such as the solar wind and the near-Earth environment (e.g., Podesta 2017; Greco et al. 2018;Fadanelli et al. 2019;Pecora et al. 2019;Gingell et al. 2020;Khabarova et al. 2021, and references therein). The characterization of current structures in turbulent plasmas is of particular interest not only because magnetic reconnection and/or different dissipation processes can occur inside (or close to) these regions, thus enabling energy conversion and plasma heating (e.g., Gosling and Phan 2013;TenBarge and Howes 2013;Osman et al. 2014;Zhdankin et al. 2014;Navarro et al. 2016;Grošelj et al. 2017;Matthaeus et al. 2020;Agudelo Rueda et al. 2021, and references therein), but also because reconnection processes occurring within such structures can in turn feed back onto turbulence itslef by playing a major role in the scale-to-scale energy transfer (e.g., Carbone et al. 1990;Cerri and Califano 2017;Loureiro and Boldyrev 2017;Franci et al. 2017;Mallet et al. 2017;Camporeale et al. 2018;Dong et al. 2018;Vech et al. 2018;Papini et al. 2019). E-mail: manuela.sisti@univ-amu.fr In order to determine the physical behavior of coherent current structures it is of foremost importance to understand how they manifest within the (turbulent) magnetic-field dynamics. This task can be separated into two main inquiries. 1) On the one hand, one must determine the current structure geometry by defining the three characteristic scale lengths of such structures, usually called thickness (the smallest), width (the medium one) and length (the biggest). Such definition of characteristic lengths was employed, for instance, by Zhdankin et al. (2013) to characterize the current structures emerging in "reduced-MHD" simulations of plasma turbulence. While in a simulation we are always in the conditions to precisely define the geometry of a current structure, this is less obvious when it comes to in-situ satellite data. A spacecraft can indeed cross a coherent current structure during its path, but it has no information about the overall geometry of it, and even a multi-spacecraft fleet can only measure the spatial variation of physical fields on scales which are in general much smaller than those of any coherent current structure. 2) On the other hand, it is also of key importance to understand how (and if) local magnetic-field configurations within the above mentioned current structures are systematically different from the magnetic configuration that belong to the rest of the (turbulent) environment. While such a characterization can be achieved by a number of procedures, from now on we will focus on the "Magnetic Configuration Analysis" (MCA) method proposed by Fadanelli et al. (2019). The MCA method consists of a modification of existing techniques that have been previously employed to investigate the local configuration of magnetic field (namely, the Magnetic Directional Derivative by Shi et al. (2005) and the Magnetic Rotational Analysis by Shen et al. (2007); see Shi et al. (2019) for a review of these different techniques). Con-Article number, page 1 of 20 arXiv:2107.14130v2 [physics.plasm-ph] 28 Oct 2021 A&A proofs: manuscript no. main trary to measures of current structure geometry, the analysis of magnetic configurations can be performed on data from plasma simulations as well as on data coming from multi-spacecraft missions. For instance, in Fadanelli et al. (2019) the authors apply the MCA technique on long intervals of data collected by the Magnetospheric Multiscale (MMS) mission (see Burch et al. (2016)) determining that it is possible to obtain statistics of local configurations developing in the outer magnetosphere, in the magnetosheath, and in the near-Earth solar wind.
In the first part of this work, we develop and describe three methods by which it is possible to estimate some aspects of a current structure geometry starting from local measurements. We test these methods on the structures emerging in 3D-3V Hybrid-Vlasov-Maxwell (HVM) simulations of plasma turbulence (three spatial dimensions plus three-dimensional velocity space, with kinetic ions and fluid electrons). By comparing the results obtained from these three local methods with those resulting from a non-local approach, we show that it is indeed possible to estimate the overall shape of a current structure by employing using only local measurements.
The second part of this work investigates the physical characteristics of current peak regions forming in three different threedimensional Hybrid-Vlasov (3D-3V HVM) simulations. In particular, we analyze the different features of magnetic configurations in the three different simulations and we make a comparison between them and those obtained analyzing the MMS observational data of plasma turbulence measured on December 9th 2016 Phan et al. (2018); Stawarz et al. (2019). Indeed, one simulation setup that we include in this work is a three-dimensional equivalent of the 2D setup employed by Califano et al. (2020), which proved capable of qualitatively reproducing the turbulent and reconnection regime observed during that period. This paper is organized as follows. In Sections 2 and 3 we present the main features of the HVM simulations of plasma turbulence that are employed here, including an overview of the turbulent spectra that we obtain in the different cases. In Section 4.1 we give a precise definition of "current structure" in our simulations and define two non-local/overall shape factors called planarity P and elongation E to better highlight the 3Dshape of a current structure. In the same Section we clarify what we mean by "magnetic configuration" and how MCA defines the two shape factors planarity P and elongation E to characterize any such configuration. Then, in Section 4.2 we present three different methods by which it is possible to convert purely local measures into an estimate of E and P of current structures, and show their effectiveness. In Section 5 we investigate the physical features of the current regions forming in the 3D simulations we make use here and perform, using our simulations, the analysis proposed by Fadanelli et al. (2019). By performing this analysis we can show that magnetic configurations inside the coherent current regions behave differently with respect to those in the rest of plasma. Moreover, we show how the statistical distributions of these quantities obtained by numerical simulations well reproduce the ones obtained by analyzing MMS data of December 9th 2016 Phan et al. (2018); Stawarz et al. (2019). Finally, we summarize the results obtained and their importance in Section 6.

Simulations
In this paper to investigate the emergence of (coherent) current structures and characterize their nature, we make use of kinetic numerical simulations of plasma turbulence performed with HVM model with kinetic ions and fluid neutralizing elec-trons with mass Valentini et al. (2007). The corresponding set of equations is normalized using the ion mass m i , the ion cyclotron frequency Ω ci , the Alfvén velocity v A and the ion skin depth d i = v A Ω −1 ci . As a result, the dimensionless electron skin depth is given by the electron-to-ion mass ratio, d e = √ m e /m i . The ion distribution function f i = f i (x, v, t) evolves following the Vlasov equation that in dimensionless units reads giving the number density n and the ion fluid velocity u i as moments of f i . Then, the electron fluid response is given by a generalized Ohm's law for the electric field E including the Hall, diamagnetic and electron-inertia effects: In the electron inertia term we approximate 1/n = 1, in normalized units, for the sake of computational simplicity. Furthermore, J = ∇ × B is the current density (neglecting the displacement current in the low-frequency regime). In equation (2) we assume an isothermal equation of state for the electron pressure, P e = nT 0e , with a given initial electron-to-ion temperature ratio T 0e /T 0i = 1. The initial ion distribution function is given by a Maxwellian distribution with corresponding uniform temperature. Finally, the evolution of the magnetic field B is given by the Faraday equation: We take an initial magnetic field given by a uniform background field, B 0 e z where B 0 = 1, with a superimposed smallamplitude 3D perturbation, δB = δB x e x + δB y e y + δB z e z , computed as the curl of the vector potential, δB = ∇ × δA. In particular, δA is given by a sum of sinusoidal modes with random phase in a limited wavevector interval corresponding to the largest wavelengths admitted by the numerical box. In Table 1 we list the simulation parameters of four different simulations as follows: the box spatial size, the corresponding number of grid points, the wave vector range of the initial perturbation, the root mean square amplitude of each component of the perturbed magnetic field and the ratio between the root mean square of the magnetic perturbation and the equilibrium field. We sample for all simulations the velocity space using 51 3 uniformly distributed grid point spanning [−5v th,i , 5v th,i ] in each direction, where v th,i = β i /2 v A is the initial ion thermal velocity and β i = 1. We set a reduced mass ratio m i /m e = 100.
In the following we will refer to the 3D-3V simulations (the first three listed in Table 1) using the following names: "SIM A", "SIM B-wf" (weak forcing scenario) and "SIM B-sf" (strong forcing scenario). We also make use as a reference for comparison of a 2D-3V hybrid Vlasov-Maxwell simulation, namely "SIM 2D" (last one in Table 1). The reasons behind the choice of simulation parameters are discussed in the following section, providing also a brief overview of turbulence evolution and properties. As summarized in Table 1 each simulation initially injects fluctuations' energy in a different wavenumber range, with different root-mean-square (rms) values of the initial magnetic fluctuations exploring a different amount of sub-ion-gyroradius scales. In particular, SIM B-wf and SIM B-sf inject energy only slightly above the ion characteristic scales, and are able to excite a wide range of kinetic scales before reaching their dissipation scale (i.e., resolving about a decade of clean sub-ion range turbulence before entering the dissipation-dominated scales). In fact, as done in Califano et al. (2020) for the 2D-3V case, the aim of both SIM B-wf and SIM-sf is to mimic the conditions possibly encountered in the Earth's magnetosheath, past the bow shock, where the occurrence of electron-only reconnection has been observed Phan et al. (2018); Stawarz et al. (2019). For completeness, we include a simulation (SIM A) where energy is injected at larger scales including part of the final MHD turbulent range and a simplified two-dimensional one (SIM 2D) covering an even larger range.
We recognize three phases in the temporal evolution of all simulated systems: (1) an initial phase, (2) a transition phase and (3) a fully developed turbulence. At the beginning (1) the energy injected at large scales starts to cascade towards smaller and smaller scales. In the transition phase (2) coherent structures, i.e. regions where the current density peaks (it usually happens at about one eddy-turnover time) start to form. At the end, a fully developed turbulent state (3) is reached, with a complex nonlinear dynamics with coherent structures continuously forming, merging and/or being destroyed. We will indicate the times at which the 3D simulations reach saturation (maximum of the root mean square of the current density) as t sat A (SIM A), t sat B−w f (SIM B-wf) and t sat B−s f (SIM B-sf); at these times turbulence is fully-developed (roughly speaking saturation times correspond to about 3 eddy-turnover times). From now on, unless explicitly stated otherwise, all the analysis will be carried out at saturation times.
In Figure 1 we show the reduced one-dimensional, kaveraged, magnetic energy spectra versus k ⊥ for the simulations listed in Table 1. The spectra, for visual purposes, have been normalized so that they overlap at k ⊥ d i ≈ 2. The power laws k −5/3 ⊥ and k −8/3 ⊥ are displayed as references for the MHD and kinetic range, respectively. At "fluid" perpendicular scales, in the range 0.1 ≤ k ⊥ d i ≤ 2, a power law close to (but slightly steeper than) −5/3 is visible in both SIM A and SIM 2D. The spectral slope evaluated at k ⊥ d i 1 for these two simulations (excluding the first modes, which are affected by the initial condition), is indeed in the range −1.9 α −1.8 (the actual value slightly depends on the exact wavenumber range on which we perform the fit). On the other hand, SIM B-wf and SIM B-sf do not include enough MHD scales to draw any conclusion regarding the emerging spectral slope for k ⊥ d i < 2. While a difference between the observed slope in SIM 2D and the expected −5/3 power law may be related to the reduced dimension that pre-vents to properly establish critical balance Schekochihin et al. (2009), this discrepancy is also observed in the spectrum emerging from the 3D case (SIM A, whose magnetic-field spectrum seems to agree very well with the corresponding spectrum of SIM 2D). As noticed also in previous 3D-3V simulations Cerri et al. (2019), whether this feature in SIM A is due to the limited extent of the MHD range or not is unclear, and will have to await further investigation by means of even larger 3D simulations. At "kinetic" scales in the range k ⊥ d i > 2, a power law fairly consistent with −8/3 clearly emerges only in SIM B-sf (a hint of such power law is visible also in SIM A, but the limited extent of the sub-ion range does not allow it to develop over a broad range of scales before the numerical dissipation takes over). Such a slope would be consistent with an "intermittency corrected" kinetic-Alfvén-wave cascade occurring at sub-ion scales Boldyrev and Perez (2012), and previously reported via 3D-3V hybrid-Vlasov simulations Cerri et al. (2017bCerri et al. ( , 2018. The other two simulations, SIM 2D and SIM B-wf, exhibit a power law steeper than −8/3. A magnetic-field spectrum close to k −3 ⊥ in the kinetic range was already observed in previous 2D hybrid-kinetic simulations (see e.g., Ref. Cerri et al. (2017a) and references therein), sometimes attributed to the effect of a reconnection-mediated twodimensional cascade at sub-ion scales Cerri and Califano (2017); Franci et al. (2017) (which, however, may continue to hold also in three spatial dimensions Loureiro and Boldyrev (2017); Mallet et al. (2017), although a definitive evidence of such regime has been elusive in 3D kinetic simulations performed thus far Cerri et al. (2019); see also Agudelo Rueda et al. (2021) for a more recent attempt). For what concerns SIM B-wf, on the other hand, a steeper spectrum may be attributed to a weaker cascade associated to the lower amplitude (with respect to SIM B-sf) of the injected fluctuations, rather than to dissipation effects associated to ion-heating and/or fluctuations' damping mechanisms in the β ≈ 1 regime considered here (see e.g., Refs. Told et al. (2015); Sulem et al. (2016); Arzamasskiy et al. (2019) and references therein). In fact, in SIM B-wf there is no evidence of a clear cascade developing along the direction parallel to B 0 (not shown here). However, a detailed analysis of fluctuations' spectral properties and associated ion-heating mechanisms (and how they change between different simulations) is beyond the scope of this work and will be reported elsewhere.
In Figure 2 we show for a) SIM A, b) SIM B-st, the 3D rendering of the current density magnitude in shaded isocontours (we remind the reader that SIM B-st employs a simulation box size that, compared to the domain of SIM A, is about one third in each direction). In blue the regions of grid points which overcome a specific threshold, defined in the next Section, on the current density J th . Fig. 1: For our four simulations one-dimensional k -averaged magnetic energy spectra versus k ⊥ , normalized so that they overlap at k ⊥ d i ≈ 2.

Definitions
In order to analyze the physics of current structures, we first define the procedure aimed at identifying a current structure as well as the parameters to be adopted in order to characterize their shape. Alongside with such definitions, we will also discuss the concept of "magnetic configuration" and the meaning of its shape, a concept which will be thoroughly applied in the following sections.
Similarly to what is done in Uritsky et al. (2010); Zhdankin et al. (2013), we adopt the formal definition of "current structure" as any connected region where the current-density magnitude, J ≡ |J|, is above a threshold value for all internal points, and the boundaries do not touch any other such zone. In particular, we stress that the specific definition of such threshold may be somewhat arbitrary. In this work, consistently with Uritsky et al. (2010) and Zhdankin et al. (2013), we choose a current threshold defined by J th = J 2 + 3σ(J 2 ), where σ(J 2 ) = J 4 − ( J 2 ) 2 and ... denotes spatial average over the whole simulation box. We want to remark that, by this procedure, any isolated point with above-threshold current density is not recognized as a current structure by its own, but it can be identified as part of some other current structure if the two are within a minimum distance of one grid-point from each other. For shorthand, we call "center" of a structure the point within the structure domain for which J is maximum.
To characterize the spatial extension (or "geometry") of a current structure, we introduce three characteristic lengths max , med and min . In particular, the "length" max indicates the longest among the structure's dimensions (this goes to infinity in 2D systems), the "thickness" min is the smallest dimension, and we call "width" med the intermediate characteristic length.
To determine these quantities, we first compute the eigenvectors of the Hessian matrix of the current density at the center of each structure. Then, we consider the plane perpendicular to the direction of least variation (which is the one associated to max ). In this particular plane we have a direction of strongest variation for the current density (associated to min ) and a direction of weakest variation (associated to med ). We define the thickness min as the maximum distance between two points with J ≥ J th along the direction of strongest variation. On the contrary the width med is the maximum distance between two points with J ≥ J th on the plane. Finally, the length max is defined as the maximal distance between two points belonging to the same structure without restricting to special planes. We note here that the definition of the thickness min is different from that of Zhdankin et al. (2013) where the use of the full width at half maximum of the interpolated profile of J is preferred. Here we decide to define the thickness using the maximum distance between two points with J ≥ J th along the direction of strongest variation, in order to remain coherent with the definition of width and length. In the 2D case we have only thickness and width, computed as above. In Appendix A we clarify what is the impact of slightly different definitions of min on our analysis.
Once that the "geometry" has been clarified, it is immediate to define also the "shape" of current structures. In a current structure, we call "shape" the quantity determined by two "structure shape factors" defined as: planarity (current str.) These two parameters measure the tendency of the current structure to squeeze toward some elongated form (E → 1) or to flatten (P → 1). In this way, we can represent each structure's shape in the plane E − P, and eventually classify shapes into categories. In particular, we adopt the following nomenclature for certain characteristic shapes: (i) "pseudo-spheres" (low E, low P), (ii) "cigars" (high E, low P), (iii) "pancakes" (low E, high P), (iv) "knife blades" (high E, high P). In Figure 3 we show a schematic representation of a current structure with its characteristics lengths: min , med and max . We note that the three quantities of thickness, width and length, here presented, are enough to define the shape of compact structures in 3D. In case the holes would appear in current structures, new parameters should be added as a filling fraction or the hole dimensions. Such upgrade will be investigated in future work. By "magnetic configuration" we mean the local characterization of the B field that can be inferred at any point in our system, even far from a current structure or its peak, by inspection of the (symmetric) tensor N ≡ [∇B] · [∇B] T /B 2 Fadanelli et al. (2019). Because of its symmetric nature, N can be defined in terms of its orthogonal eigenvectors and corresponding eigenvalues: we denote these eigenvalues by σ max , σ med , σ min . The quantities 1/ √ σ min , 1/ √ σ med , 1/ √ σ max can be understood as proportional to the three characteristic variation lengths. Here and elsewhere, when considering the eigenvalues, it is important to stress that they are only proportional to the characteristic variation lengths, not directly comparable, due to the presence of normalization factors. What is really significant is the ratio between eigenvalues, as we will see in the next. Indeed, taking the ratio between eigenvalues the normalization factors cancel each other out.
As with current structures, a "shape" can be defined for magnetic configurations. The idea is to deal with 1/ √ σ min , 1/ √ σ med and 1/ √ σ max in the exact same way we did with characteristic lengths, and define two "configuration shape factors" as follows: We note here that the "configuration shape factors" are local quantities, which can be computed at any spatial location, including points owning to a current structure. On the contrary the "structure shape factors" provide an overall/non-local picture of a current structure. As E and P do for current structures, E and P do measure the tendency of magnetic configurations toward elongation (if E → 1) or squashing into sheets (if P → 1). Similarly, the position of a configuration in the E -P plane allows to recognize the local shape of B in one of the following four categories: (i) "pseudo-sphere", (ii) "cigar", (iii) "pancake" or (iv) "knife-blade". We remind the reader that the work by Fadanelli et al. (2019) is focused only on satellite data and thus on local measurements, providing point by point value for E and P along spacecrafts' trajectories. On the contrary here, taking advantage of the 3D data of a simulation, we can also give an overall/non-local picture of current structures (via E and P).

The shape of current structures
The first goal of this Paper is to determine whether it is possible to estimate the "non-local" overall shape of a current structure by local measurements, i.e. if there is some procedure to convert purely local measurements of certain quantities into a reliable estimate of the "actual" elongation and planarity of a given structure. We point out here that we consider the estimate of E and P as given by our "non-local" overall method in equations (4) and (5) to be used as "reference values" for the estimates obtained via the following local methods: -"HJ" method. This approach is based on the Hessian matrix of J calculated at the center of each current structure. If we suppose that ratios of the eigenvalues of this matrix, h max , h med and h min , can reproduce the ratios between max , med and min , then the shape factors of the structure can be estimated as: where the subscript "cen" indicates that the values are calculated at the structure's center. This method, we note, closely resembles the one adopted by Servidio et al. (2009). -"NB" method. This approach aims at inferring the "structure shape factors", E and P, from the "magnetic-configuration shape factors", E and P, evaluated through the matrix N at the center of a given current structure: where the subscript "cen" indicates that the values are calculated at the structure's center only, and σ max , σ med , σ min are eigenvalues of the tensor N which evaluates the magnetic configuration. it is worth noticing that, while the N tensor is well defined everywhere, for the sake of comparison with the other methods we consider here only the value at the center assuming it as representative for the whole structure. However, to easily make a comparison with satellites data and/or to increase the dataset analyzed, one can apply this method even to points different from the central one. Since magnetic field and current density are strongly related (in our model, the current density is the curl of the magnetic field), we also expect E NB and P NB to resemble E HJ and P HJ . -"AV" method. This approach considers an average of the magnetic configurations occurring inside each structures. That is, shape parameters are estimated as: where ... str is the average over all points belonging to a single current structure. This procedure, we note, resembles the NB method but can be carried on without knowing where the center of a current structure is located.

Method comparison
Our goal here is to verify the accuracy of the different methods HJ, NB and AV in estimating E and P, that are obtained looking at the overall/non-local shape of the current structure and are thus considered our reference values. To establish how accurate the different methods detailed in Section 4.2 are, in this Section we apply them to SIM A, which employs the largest physical domain (and thus is the simulation developing the largest number of current structures to test). To this end, we compute elongation and planarity of every current structure, then compare them with their estimates obtained by the HJ method, NB method and AV method.
In Figure 4 we show how E and P are distributed when defined through max , med , and min (Equations 4 and 5), while in Figure 5 we report the results of the estimates from HJ, NB, and AV methods (i.e. distributions for E HJ and P HJ , E NB and P NB , E AV and P AV , respectively, alongside with their ratios to the reference "non-local" overall E and P values for each structure). The physical picture that emerges from Figure 4 is consistent with highly elongated current structures. Most of the occurrences are concentrated at elongation in the interval 0.8 − 0.95 and planarity between 0.5 and 0.8 (in particular, for mean, median and standard deviation see first two rows of Table 2) suggesting "knife blade" 3D configurations. The results from HJ and AV methods are in very good agreement with the overall occurrence distribution of E and P, as shown in the top-row panels of Figure 5. For what concerns method NB, it gives the same physical picture of the other local/non-local methods, with highly elongated and mostly planar current structures (this is particularly clear from panel (e)), but from panel (b) we also note a tendency to return higher values of planarity and elongation (indeed the distribution is a bit shifted towards the upper-right corner). The near equivalence of the three methods in evaluating the characteristic shape of the current structures is further confirmed in the bottom row of Figure 5. More specifically, Figure 5 (d) shows that the distribution of E/E HJ is peaked around 0.9, and that P/P HJ is peaked around ∼ 1, meaning that for the majority of current structures the true elongation and planarity are almost the same of the ones obtained using HJ method (for mean, median and standard deviation see 3rd and 4th rows Table 2). The same agreement is shown in Figure 5 (e) between the control values and the NB method estimates, as shown by mean, median and standard deviation reported in Table 2 (5th and 6th rows). Finally, Figure 5 (f) shows that the AV method performs even better as compared to other local methods and this is confirmed by the mean, median and standard deviation values reported in Table 2, 7th and 8th row.
Despite the overall picture provides a very good agreement between reference and local methods, we note the presence of a slightly but persistent underestimation of the planarity using the local methods with respect to the reference one, underlined by the mean values reported in Table 2. Despite this general effect, we deduce that by using local methods, and in particular those employing the N tensor, it is possible to correctly estimate the actual shape of current structures. The main advantage of using the N tensor is that it can be applied to any point and not necessarily on the center of current structures. Even more importantly, the AV method, which is also the one which performs the best, can be applied with slight modifications to satellite data, once that a proper current threshold is defined.

Investigating physical characteristics of current structures
In this section we investigate the physical features of the structures forming in the three different 3D simulations. We will investigate the magnetic configuration shape factors (planarity and elongation), the characteristic eigenvalues (proportional to the characteristic scale lengths) and the orientation of magnetic field and current density. We will use the N tensor and its eigenvalues, but without restricting our analysis to points owning to the current structures, as for NB and AV methods, thus we will essen- Fig. 4: Occurrence distributions of elongation E versus planarity P (as computed using the "non-local" overall method discussed in Section 4.1) for SIM A.  tially apply the MCA method as described first in Fadanelli et al. (2019). The aim is to determine what are the statistical characteristics of magnetic configurations that emerge in our simulations and whether there is any appreciable difference between the configurations found inside and outside current structures. Thus, we will consider separately 1) all points belonging to current structures, 2) the same number of points as in the first case, but belonging to an uniform sampling of the simulation box. These last points are called "generic". This type of analysis further expands the methodology employed by Fadanelli et al. (2019) on satellite data, where, however, the statistical properties of magnetic configurations have been investigated without distinguishing current structures from the rest of the plasma. In particular, the analysis of "generic points" can be considered the equivalent of a continuous sampling of satellite data, as done by Fadanelli et al. (2019). We expect that the uniformly picked set would reproduce data collected by satellite along 1D trajectories. Indeed, in the uniform sampling we picked one point every three along each direction. Because the typical dimensions of a current structure are in average big enough to include more than three grid points in at least two directions, for each structure several points are collected, as it is the case for a continuous satellite sampling. Investigating the possible difference between an uniform sampling and synthetic 1D satellite trajectories is out of the scope of the present manuscript. Furthermore, we compare the above analysis performed on our 3D-3V simulations with an analogous analysis applied to two intervals of high-resolution (burst) data which have been collected as the MMS spacecraft encountered a turbulent magnetosheath region just downstream the bow shock. These same intervals were previously considered by Stawarz et al. (2019) for a complete analysis of turbulence. In order to be directly comparable with our simulations, a three-step selection on the abovementioned MMS data has been applied so that the magnetic configuration analysis has been performed only on those data points for which (i) the computation of N is precise enough that at least two eigenvalues are well determined, (ii) β ∈ [0.3, 3] and (iii) the resolution attained by MMS data is comparable with that of the numerical simulations (after these selections, roughly 20% of original data are kept; for details, see Appendix B).

Analysis of shape factors for magnetic configurations
In Figure 6 we show the occurrence distribution of E vs P in our set of 3D simulations, namely SIM A, SIM B-wf, and SIM B-sf (left, center, and right column, respectively): in the top row for generic simulation points (which can be considered as the equivalent of the continuous sampling in the methodology developed by Fadanelli et al. (2019)), in the bottom row only for points be-longing to what is identified as a "current structure" (i.e., those regions where the current density is above J th ).
Occurrence distributions which consider only points belonging to current structures (bottom row) or considering generic points (top row) look very different. In particular, if we consider generic points we get a distribution spread in planarity between ∼ 0.2 and ∼ 0.9 with a peak around ∼ 0.6 − 0.7. Instead, if we consider only points which belong to current structures we get an increase in planarity with the distribution peaking around ∼ 0.8 − 0.9 in the y-axes. The fact that the distributions of magnetic configurations obtained via generic simulation points does not provide a reliable estimate is even more evident when comparing the results for SIM A presented in Figure 6 with those obtained in Section 4.3 for the same simulation (i.e., see Figures 4 and 5(a-c)). We further point out that for SIM Bwf this discrepancy between the two distributions (viz., generic points vs only points belonging to current over-densities) seems to be particularly emphasized. In the following we will give the mean, median and standard deviation for the corresponding 1Ddistributions of planarity and elongation.
The statistical features of the magnetic configurations can be appreciated in Figure 7 where we show the 1D normalized occurrence distribution of E and P for the three simulations SIM A, SIM B-wf and SIM B-sf for generic points (magenta line) and for points belonging to current structures (black line). We have superposed to these distributions obtained from simulations, the ones obtained using MMS satellites data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019) (see Appendix B for details).
As in the color-maps of Figure 6, the histograms in Figure 7 shows an appreciable statistical difference in planarity between generic points (magenta lines) and those located inside current structures (black lines). Indeed, if we consider generic points we obtain a P distributions less skewed towards P ≈ 1, with a peak around 0.6 − 0.8 (please, note the log-scale for the y-axis). On the contrary, the normalized histogram of P for points owning to current structures peaks around 0.8 − 1.0. For all three simulations the values of mean, median and standard deviation are reported in Table 3 (1st row for generic points and 2nd row for current structures) and turn out to be the same. For the elongation, the behavior changes less significantly when considering generic points or current structures. Indeed, the peaks of the occurrence distributions are always close to 1. Also in this case we report mean, median and standard deviation in Table 3 (3rd row for generic points and 4th row for current structures). Summarizing, the emerging picture indicates a majority of "blade-like" magnetic configurations. Commenting on the different behavior between generic points and those located in current structures, this is not surprising since the points located in current structures belong to regions which tend to have a specific shape and a common behavior. In particular, we note that they have an high planarity confirming the intuitive picture of nearly 2D current sheets. On the contrary for generic points, picked up almost from everywhere, also from regions where the current density is low, the magnetic configuration has a different behavior.
Finally, a major result is the behavior of the distributions for satellite data in agreement with the results of our simulations for generic points. In particular, the agreement is the most evident for SIM B-sf. This agreement is a very good result but not surprising for the following two reasons. First, since there is no selection on the values of J for the analysis performed on MMS time series, it was expected that any agreement with such analysis would have involved the sub-set of "generic points" from our simulations. Second, simulation SIM B-sf has a setup which is   (2020), which already proved capable of qualitatively reproducing the turbulent and reconnection regime observed during that period. Thus we expected to find similar features in the configuration shape factors comparing SIM B-sf and these particular MMS intervals.

Analysis of eigenvalues' occurrence distribution and characteristic scale lengths, in simulations and satellite data
In Figure 8 we show the normalized occurrence distribution of N's eigenvalues: a) σ min , b) σ med , c) σ max for SIM A, SIM B-wf and SIM B-sf. We recall that these eigenvalues are interpreted as proportional to the squared inverse of the local variation length (see Section 4.1). We superpose to these distributions, obtained from simulations, the ones obtained using satellites data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019) (with the additional selection discussed at the beginning of this Section; see Appendix B for Fig. 7: Normalized occurrence distributions of P and E for the three simulations SIM A, SIM B-wf and SIM B-sf, distinguishing between statistics on generic points (magenta line) and points which belong to current structures (black line). We superpose to these distributions, obtained from simulations, the ones obtained using satellites data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. details). In the following, we note that the MMS eigenvalues are normalized using the local value of d i .
Let us consider first the occurrence distribution of the smallest eigenvalue (which is related to the largest characteristic length of the magnetic configuration). This is shown in column a) of Figure 8. For all three simulations there is almost no difference between the distributions for generic points (magenta line) and for points belonging to current peaks (black line) related to this eigenvalue. Indeed the two lines are almost superposed. This means that the typical largest scale length of a "generic" region does not differ significantly from that of a structure where the current peaks. Now let us consider instead the occurrence distributions for the other two eigenvalues (related to the median and shortest variation length of the local magnetic configurations; shown in column b) and c) of Figure 8, respectively). In all the simulations, the occurrence distributions of these eigenvalues Fig. 8: Normalized occurrence distribution of N's eigenvalues: a) σ min , b) σ med , c) σ max , for SIM A, SIM B-wf and SIM B-sf, respectively). In magenta, distributions for generic points, in black, for points belonging to current structures. We superposed to these distributions, obtained from simulations, the ones obtained using satellites data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. Note the logarithmic scale on the x-axis. evaluated at generic points significantly differ from the corresponding distributions that are obtained considering only those points belonging to current structures. In particular, the difference is more pronounced for the smallest length scale, i.e. the one associated to the largest eigenvalue. In all cases, when compared to those obtained using generic points, the distributions obtained using only the points belonging to current structures shift towards larger σ. This means that where the current attains values higher that J th , the typical scale length of variation of the magnetic field are smaller (as it is somewhat expected, since the current is the curl of the magnetic field). Comparing SIM B-wf and SIM B-sf, which we remind the reader are identical except for the amplitude of the injected perturbation, one finds that the peak of the occurrence distributions is located at larger values of log σ med and log σ max for simulation SIM B-sf rather than for SIM B-wf. Thus, simulation SIM B-sf develops current structures with typical scale length smaller than those of the other simulations. This is true also for generic points, even if less pronounced. Such feature is likely a consequence of the fact that the two simulations develop different turbulent regimes (see discussion in Section 3). In fact, by injecting larger-amplitude magnetic fluctuations in SIM B-sf than in SIM B-wf, a shallower sub-ion-scale spectrum develops and, consequently, a larger amount of turbulent power reaches the smallest scales (cf. the magneticfield spectrum in Figure 1). Whether this is actually due to a weak-like versus strong-like turbulent regime (i.e., if it could be a consequence of whether critical balance and/or dynamic alignment establishes or not) or to other effects (e.g., development of a significant amount of electron-only reconnection events) is beyond the scope of this work and will have to await further investigation, as well as additional simulations. For current structures, we report the mean, median, standard deviation values in Table  3, 5rd-8th rows, comparing quantitatively SIM B-wf and SIM B-sf.
Concerning the comparison with the distributions extracted from MMS data Stawarz et al. (2019), we note a very good agreement with the distributions for generic points for SIM Bst. As for the previous Section 5.1, this very good agreement with SIM B-sf was somewhat expected since a similar simulation setup in the 2D case already proved to be capable to qualitatively reproducing the turbulent and reconnection regime observed in these intervals Califano et al. (2020). Thus, we proved once again that a relatively large-amplitude injection close to the ion-kinetic scale is able to reproduce several features that are observed by MMS in the turbulent magnetosheath past the bow shock-in particular, the development of similar 3D structures. Although such setup constitutes a quite interesting hint, the actual mechanism that could be behind such peculiar injection (e.g., the bow shock itself or the subsequent occurrence of micro-instabilities) is still unclear and requires further investigations, especially from the observational point of view.

Analysis of the orientation of magnetic field and current density
We analyze the orientation of the direction e min along which we measure the smallest eigenvalue with respect to the direction of the local magnetic field and current density by computing |b·e min | and |j · e min |, respectively. Here b and j are the unit vectors of the local magnetic field and current density. As it was done in previous Sections, this calculation is performed both on generic simulation points and also restricting only to points which belong to current structures. In Figure 9 we show the normalized occurrence distribution of |b · e min | (left column) and |j · e min | (right column) for all three simulations for generic points (magenta) and for points belonging to current structures (black). We superpose to these distributions those obtained by using satellites data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019) (see Appendix B for details).
Both the local magnetic field and j are well aligned to e min for all three simulations. In particular, we note that the alignment is stronger for points which belong to current peaks, which means that for these configurations the alignment between magnetic field and e min (and the same for j) is the best one. The strong alignment between j and e min is expected since by definition the derivatives of B perpendicular to e min are the strongest ones. Instead, the alignment between the magnetic field and e min is less obvious and could depend on the specific environment that is being considered and/or on the initial parameters of our simulations. In general, we expect that the local magnetic field and the current density tend to align only when there is a significant component of "guide field" in the structure (opposed to those configurations where B vanishes within the current structure, as, for instance, in the typical setup that is employed to study magnetic reconnection without a guide field). Moreover, the good alignment between e min and B could also be due to Beltramization of the flow, namely the alignment between current and magnetic field, which is typical of small-scale turbulent structures (see e.g. De Giorgio et al. (2017)).
Also in this case, the comparison with the observative data from MMS Stawarz et al. (2019) is very good. In particular, we note that for the alignment |b · e min | the best agreement between MMS data and simulations "generic points" is found with SIM B-wf rather than with SIM B-sf. This is probably due to the fact that |b · e min | is strongly affected by the presence of a guide field within the current structure. In fact, if δB/B 0 is small enough as it is for SIM B-wf, the guide field is less distorted by the turbulent dynamics, and also the direction of weakest variation of the emerging current structures will align better with b. On the other hand, when large δB/B 0 fluctuations are injected, there is a significant distortion of the background magnetic field. As a result, the emerging structures can be embedded in magnetic shears where, with respect to the plane perpendicular to e min , there is a weak (or vanishing) guide field.
We note, as in the discussion of Figures 7 and 8, that the agreement of the MMS distribution with the "generic points" rather than with the dataset of points belonging to current structures was to be expected. Indeed, there is no selection on the values of J for the analysis performed on MMS time series.

Generic points or "background"
We briefly explain what implies to perform a statistical analysis on "generic" points of simulations or, equivalently, on long continuous sampling in observation data, referring to the methodology proposed in Fadanelli et al. (2019). First of all, we want to make sure that the ensemble of "generic" points picked in our analysis constitute a somewhat statistically representative subgroup of the whole grid points of the simulation, while simultaneously being comparable to the total number of points that belong to the current structures in such simulation (which is instead the sub-set used when restricting the analysis to these structures). Let us consider, for instance, SIM A, in which the total number of grid points is 352 * 352 * 198 ∼ 2.5 × 10 7 . In this case, the sub-set of "generic" points has been created by considering one point every 27 (corresponding to a collection of ∼ 10 6 points), i.e. around ∼ 4% of the original set. Analogously, the total number of points which overcome the threshold on current density at time t sat A are ∼ 1.1 × 10 6 (i.e., roughly 4% of the total). Moreover, in our sub-set of "generic" points those which overcome the threshold on current density are 4.8 × 10 4 , again roughly 4% of the sub-set. Thus, the sub-set of "generic" points that is being considered should adequately represent the ensemble of points of our simulation box. As we just discussed, the points which constitute the current structures are always a small percentage of the total number of points in a set or sub-set (i.e., the filling factor of the current structures in a volume is usually very small). Therefore, a natural question to ask is the following: can the behavior of a such small sub-group emerge when we consider histograms and plots which refer to "generic" points? (or, equivalently, when we consider a long continuous sample in observation data?) Based on our analysis, the answer to this question is no. This can be seen in Figure 10, which shows the occurrence distribution of elongation E versus planarity P for a) only those "generic" points that do not belong to current structures (roughly 4%), b) the whole sub-set of "generic" points (i.e., including those belonging to current structures). In fact, there is no noticeable difference between the two distributions. In partic- Fig. 9: normalized occurrence distribution of |b · e min | (left column) and |j · e min | (right column) for all three simulations, in magenta the distribution for generic points, in black that for points belonging to current structures. The distribution refers to times at fullydeveloped turbulence, in particular t sat A for SIM A, t sat B −w f for SIM B-wf and t sat B −s f for SIM B-sf. The distribution refers to times at fully-developed turbulence, in particular t sat A for SIM A, t sat B −w f for SIM B-wf and t sat B −s f for SIM B-sf. We superposed to these distributions, obtained from simulations, the ones obtained using satellites data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. ular, for both cases the elongation and planarity have the same mean, median and standard deviation (see Table 3, rows 9 to 12).
The behavior of the current structures cannot emerge if we consider "generic" points, i.e. the equivalent of a long continuous sample in satellite data without introducing any additional selection on the turbulent time series. Therefore, in order to sys-tematically study the statistical behavior of the shape of current structures in satellites data, we suggest to fix a threshold on the current density and to perform the analysis only on those regions that overcome this threshold. This procedure is similar to what has been done when applying the AV method discussed in Section 4.2 to our numerical simulations. A similar kind of sugges-

Conclusions
In this work we have conducted a wide-spanning analysis of overall/non-local current structures and local magnetic configurations emerging in plasma turbulence. The analysis is based on three 3D-3V HVM simulations with different box sizes, resolution and energy injection scales. We focused first on the characterization of the shape of current structures. Our "reference method" defines two parameters, elongation E and planarity P, which can be calculated in a simulation from the three characteristic dimensions of any structure ( max , med and min ). We have shown how it is possible to reliably estimate the shape of a current structure also by employing three different local methods, namely HJ, NB and AV (see Section 4.2 for details on the methods). The rigorous computation of E and P via the "reference" non-local method can be performed only on simulation data. The HJ and NB methods too can be performed only on simulation data, since they require to be applied on the center of a current structure (whose precise position can be known only on simulation data). These two "local" methods are faster than the other two (i.e., the "non-local" and the AV methods), since they require to be computed on a restricted number of points (namely, only on the local maxima of the current density). Finally, the AV method, the one in best agreement with the results of the "reference" method, can be applied with minor modifications also on spacecraft data, once a threshold on the current density has been properly defined on the time series under consideration. Based on these four methods (i.e., "reference/non-local", HJ, NB and AV), we have analyzed the distribution of shape factors (i.e., planarity and elongation) for the emerging current structures, with the result that all methods coherently find that they are composed of mainly "knife-blade"-like structures. This picture is different from the one that emerges in Meyrand and Galtier (2013) where they claim the presence of mostly cigarlike structures (they use the expression "filament like") through visual inspection of their 3D Hall-MHD simulations of turbulence. This suggests that ion-kinetic effects (i.e., beyond just the Hall term) and/or electron-inertia terms could significantly affect (and likely be required to correctly describe) the development of current structure in plasma turbulence across the transition range and at sub-ion scales.
Additionally, we studied the local magnetic configurations by performing an analysis on the simulation data similar to the one proposed by Fadanelli et al. (2019) for satellite data. In particular, we have investigated the magnetic configurations by analyzing the distribution of their planarity and elongation, of their three characteristic scale lengths, and of their orientation with respect to the magnetic-field and current-density directions. Such analysis has been performed both on all points belonging to current structures as well as for "generic" points belonging to a uniform sampling of the simulation box (the aim of this latter set being to mimic long and continuous time series from satellite data). In general, we found different results when we apply the analysis only to those points belonging to current structures or to "generic" points in the simulation domain (i.e., including, but not limited to, structures). In particular, the main statistical difference between generic points and those located inside current structure is in the results obtained for the distribution of planarity P: "knife-blade" shapes are more likely present when considering only those points where the current density is above a certain threshold J th , while the abundance of "thicker" sheets (or "ellipsoids") is enhanced when points below J th are included. The behavior of the distributions for generic points for planarity P and elongation E is coherent with the distribution we obtain if we apply the same kind of analysis to high-resolution (burst) MMS data collected from the two turbulence crossings in the magnetosheath and analyzed previously by Stawarz et al. (2019), and selected in order to fit simulations' characteristics (see Appendix B for details).
In the analysis of variation scale lengths, we found a sensible difference in the distributions for the largest eigenvalue between generic points and points within structures, suggesting that the smallest characteristic length-scale min is significantly shorter for current structures. We noted a difference in this context also between SIM B-wf and SIM B-sf, for which all the three characteristic length scales are shorter in the strong-forcing scenario (SIM B-sf) than they are when lower-amplitude fluctuations are injected (SIM B-wf). We remind the reader that in both simulations SIM B-wf and SIM B-sf the energy is injected only slightly above the ion characteristic scales, as done in Califano et al. (2020), that was able to reproduce, in a simplified 2D-3V configuration, the electron-only turbulent regime observed by MMS past the bow shock Phan et al. (2018); Stawarz et al. (2019). Concerning the comparison with the distributions extracted from these precise MMS data Phan et al. (2018);Stawarz et al. (2019), we found a very good agreement only with the distributions for generic points from SIM B-sf, thus suggesting that similar, small-scale current structures require a high level of forcing acting close to the ion-kinetic scales.
Finally, from the analysis of the local orientation of the magnetic field and the current density with the minimum variance direction e min , we found strong alignments for both fields, more pronounced for points which belong to current structures. Also in this case, the comparison with the observative data from MMS Phan et al. (2018); Stawarz et al. (2019) is very good. In particular, we note a good agreement with the distributions for generic points in SIM B-wf.
All the analysis presented in this work clearly highlights how magnetic configurations inside current structure exhibit peculiar features that can be retrieved only by considering solely those points where J attains values above a certain threshold. Indeed, we have shown in Section 5.4 that the behavior of the current structures cannot emerge if we consider "generic" points, i.e., analogous to consider a long continuous time series in satellite data without any further selection based on the values of J, since for such a sample the number of points belonging to current structures will be only a small percentage of the total. Therefore, in order to systematically study the shape of current structures in satellites data, we suggest to fix a proper threshold on the current density and consequently consider only regions that overcome this threshold for the subsequent analysis. Even better, one could apply the AV method that we have described in Section 4.2 in order to isolate different structures (after applying straightforward modifications in order to adapt the method to simple 1D time series rather than to complex 3D spatial domains).
In conclusion, the results reported in this paper would be useful not only for the analysis of turbulence simulations, but also for observative studies. Indeed, 1) we have validated the possibility to apply local methods, which are the only ones applicable on satellite data, to infer the overall/non-local shape of current structures; 2) we conjecture that imposing a proper threshold on the current density would benefit for the statistically study of current structures in satellite data; 3) we have provided an overview of local magnetic configurations emerging in different turbulence regimes, also stressing the different behavior that is found when considering exclusively points within current structures with respect to what emerges from the points belonging to the rest of the turbulent environment; 4) we have shown via such magnetic configuration analysis that, when there is a mechanism that inject relatively high-level fluctuations close to the ion-kinetic scales, our 3D-3V simulation can reproduce the structures that emerge in MMS data for the periods studied by Phan et al. (2018); Stawarz et al. (2019).

Fig.
A.1: scatter plot thickness vs current density for SIM 2D at t eddy−turnover 2D , to show the differences between the two methods in estimating thickness. The current density is computed at the center of each current structure. With "mask method" we mean the method to find thickness explained in Section 4.1, which uses the maximum distance between points satisfying J ≥ J th . On the contrary, with "FWHM method" we refer to the alternative method explained in Appendix A which computes thickness using full width at half maximum of the interpolated profile of J. In the sub-panel we show a schematic representation to explain the difference in estimations between the "FWHM method" and the "mask method".

Appendix A: Alternative definitions for structure's dimensions
In our work, we defined the thickness min as the maximum distance between two points with J ≥ J th along the direction of strongest variation of the current density, as given by the eigenvectors of the Hessian matrix. We present here an alternative definition for thickness, as proposed by Zhdankin et al. (2013), and we show how some of our results change when using this different definition. In particular, we can alternatively define thickness considering the interpolated profile of J along the direction of strongest variation given by the Hessian matrix and computing it as the width at half maximum of this profile. The two methods can produce different estimations for thickness depending on the value of the current density peak. In Figure A.1 we show the scatter plot thickness versus current density for SIM 2D at t eddy−turnover 2D : in red the estimations given by the method used in the main text (maximum distance between two points with J ≥ J th , for brevity, called "mask method"), while in blue the ones given by the method used in this appendix (full width at half maximum of the current density profile, called "FWHM method"). We chose the 2D simulations for this comparison and the eddy-turnover time rather than the saturation time just because in this simulation, at this time, current structures are not so many to make the scatter plot difficult to be interpreted. The current density is computed at the center of each current structure considered. We can see that the two estimations turn out to be different, especially at low values of current density. In the sub-panel, we show a schematic representation to explain the difference in estimations between the "FWHM method" and the "mask method". In particular, when the current density at the center of the specific structure is greater than 2J th , thickness estimated using the "FWHM" method is generally bigger than the one estimated using "mask method", on the contrary when the current density is lesser than 2J th the thickness obtained with the "mask method" is the biggest one. The two values are almost equal when J ∼ 2J th . If we use this alternative definition of thickness for computing the structure's shape factor planarity, some features change (we note that the estimation of elongation does not change since we haven't changed the definitions of width and length). In particular, in Figure A.2 we show the new distribution for E and P which emerges, which is a bit different from the one of Figure 4. In particular, to quantify, now we have: for the elongation mean∼ 0.8, median ∼ 0.9, standard deviation ∼ 0.2, instead for planarity mean ∼ 0.4, median ∼ 0.4 and standard deviation ∼ 0.2. The distribution is lesser in agreement with the ones for HJ, NB, and AV method showed in Figure 5, and the most of the structures have a planarity in between 0.2-0.6, thus are more filament-like.
The width (the intermediate length) is defined as the maximum distance between points having J > J th , in the plane perpendicular to e min and passing from the current peak of the structure. For the sake of coherence, we prefer the "mask method" for defining the thickness. In such a way, a cylindrical structure would have, as it is, coherent values for both width and thickness, independently from the intensity of the current.

Appendix B: Selection and processing of MMS data
In a number of passages in the main text, we referred to the two burst-resolution MMS data intervals analyzed in Stawarz et al. (2019), which have been chosen here for comparison with the results of numerical experiments detailed in Section 2. Choice of these MMS data intervals is motivated by the fact that our simulation set-up (in particular for SIM B-wt and SIM B-st) is similar to that used in Califano et al. (2020), which was able to reproduce the turbulent and reconnection regime observed during that period.
Throughout our analyses, MMS data relative to the magnetic field has been taken from the FluxGate Magnetometer (FGMsee Russell et al. (2015)), while the Fast Plasma Investigation (FPI -Pollock et al. (2016)) has provided the density and pressure measures which contributed to determine the ion plasma beta and ion inertial lengths. All these aforementioned fields have been interpolated onto the MMS2 magnetic field data, and averaged over the four spacecraft fleet. The J value has been obtained by applying the curlometer technique by Dunlop et al. (1988) and the N tensor was derived as in Fadanelli et al. (2019) i.e. by performing a linear estimation of the magnetic field gradient, then combining the resulting values with the four-spacecraft averaged magnetic field data.
In order to compare results from MCA applied over simulation and MMS data, we need apply three different levels of selection to the latter. In particular: 1) we select points for which at least two of the eigenvalues of the N matrix are well-determined by MMS measurements; 2) we select data with β ∈ [0.3, 3]; 3) we select data with a resolution comparable with the one of our simulations.
More in detail, the first selection has been performed exactly like in Fadanelli et al. (2019) i.e. by calculating the average interspacecraft distance S C at each instant and then setting at ( δB S C B ) 2 a minimal resolution threshold for the eigenvalues of N, being understood that any eigenvalue below such a threshold is not well-resolved. Requiring that at least two of the eigenvalues are well-resolved signifies that, on one hand, we accept uncertainty on the least eigenvalue and elongation measures but, on the other hand, we are still able to determine correctly the direction of minimum variation (this is because the corresponding eigenvector is, by construction, perpendicular to the other two, which are well determined). The requirement that two eigenvalues are well resolved is generally a good compromise between the need for precision and the difficulty of determining correctly the smallest eigenvalue, which is generally extremely small and therefore easily falls below the MMS resolution threshold. Choosing the two-well-determined-eigenvalues selection criterion for the two turbulence intervals we consider here implies that over the 99% of original data is retained by this first procedure. The second selection just described is performed over the ion plasma beta obtained by the ratio of spacecraft-averaged kinetic and magnetic pressures at each datapoint. For the turbulence intervals we consider in this work, this is the procedure which leads to the largest reduction in the availabel dataset, which gets reduced to about 20% of its original size through it.
With the third and last filtering of MMS data we intend to eliminate all those points for which simulation and satellite data have different resolutions. To compare precision of MMS data and of simulations, we introduce two quantities that we call "resolution factors" and are defined as follows: for the MMS (satellite) measurements, the resolution factor is the inter-spacecraft separation divided by ion inertial length (we remember that the only possible derivative calculation in this case follows from a linear interpolation of satellite measures) for the simulations, the resolution factor is max( 3dx d i , 3dy d i , 3dz d i ) and by this definition we intend to acknowledge that the effective minimal distance over which a HVM numerical simulation can well represent plasma physics should be several times (here we have chosen three) the distance between neighboring data points.
When considering MMS data analyzed in Stawarz et al. (2019) the resolution factor is almost always below 0.3. This resolution factor is thus comparable with the one we obtain for our simulations, in particular for SIM B-wf and SIM B-sf, for which it is 0.22, in all directions (instead, we obtain 0.9 for simulation run SIM A, corresponding to the z-direction which is the less resolved, while for the x and y-direction we obtain 0.5). Given these values, we have decided to accept all the previously selected data into a MCA procedure which is to be compared to that performed over HVM simulation results.
By the whole procedure just detailed, a dataset containing about 200 thousand points has been selected. Given these data, it is possible to obtain a valid statistic only for what we called "generic" points, since any further selection aiming to retain only samples retrieved inside current structures would leave us with no more than few thousand points, which are not sufficient for statistics in our case.