Issue 
A&A
Volume 655, November 2021



Article Number  A107  
Number of page(s)  18  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202141902  
Published online  29 November 2021 
Characterizing current structures in 3D hybridkinetic simulations of plasma turbulence
^{1}
AixMarseille University, CNRS, PIIM UMR 7345, Marseille, France
email: manuela.sisti@univamu.fr
^{2}
Dipartimento di Fisica, Università di Pisa, Pisa, Italy
^{3}
Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA
Received:
29
July
2021
Accepted:
24
September
2021
Context. In space and astrophysical plasmas, turbulence leads to the development of coherent structures characterized by a strong current density and important magnetic shears.
Aims. Using hybridkinetic simulations of turbulence (3D with different energy injection scales), we investigate the development of these coherent structures and characterize their shape.
Methods. First, we present different methods to estimate the overall shape of the 3D structure using local measurements, foreseeing an 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 (Magnetospheric MultiScale mission) dataset having similar plasma parameters.
Results. Thanks to our analysis, (1) we validate the possibility of studying the overall shape of 3D structures using local methods, (2) we provide an overview of a local magnetic configuration emerging in different turbulent regimes, (3) we show that our 3D3V simulations can reproduce the structures that emerge in MMS data for the periods considered.
Key words: methods: data analysis / methods: statistical / methods: numerical / turbulence
© M. Sisti et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Turbulent, magnetized plasmas permeate a wide range of space and astrophysical environments, and plasma turbulence naturally develops coherent structures characterized by a 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 (Zhdankin et al. 2013, 2017; Passot et al. 2014; Navarro et al. 2016; Cerri et al. 2019; Comisso & Sironi 2019, and references therein), as well as routinely observed via in situ measurements in space plasmas such as the solar wind and the nearEarth environment (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 (Gosling & Phan 2013; TenBarge & 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 itself by playing a major role in the scaletoscale energy transfer (Carbone et al. 1990; Cerri & Califano 2017; Loureiro & 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).
In order to determine the physical behavior of coherent current structures, it is of foremost importance to understand how they manifest within the (turbulent) magneticfield 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 a definition of characteristic lengths was employed, for instance, by Zhdankin et al. (2013) to characterize the current structures emerging in “reducedMHD” (magnetohydrodynamics) simulations of plasma turbulence. While in a simulation we can always 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 multispacecraft 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 magneticfield configurations within the abovementioned current structures are systematically different from the magnetic configuration that belongs to the rest of the (turbulent) environment. While such a characterization can be achieved by a number of procedures, from now on we focus on the magnetic configuration analysis (MCA) method proposed by Fadanelli et al. (2019). The MCA method consists of modifying existing techniques that have been previously employed to investigate the local configuration of the 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). Contrary 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 multispacecraft missions. For instance, in Fadanelli et al. (2019), the authors applied 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 nearEarth 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 3D3V HybridVlasovMaxwell (HVM) simulations of plasma turbulence (three spatial dimensions plus threedimensional velocity space, with kinetic ions and fluid electrons). By comparing the results obtained from these three local methods with those resulting from a nonlocal approach, we show that it is indeed possible to estimate the overall shape of a current structure by employing only local measurements.
The second part of this work investigates the physical characteristics of current peak regions forming in three different threedimensional HybridVlasov (3D3V 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 by analyzing the MMS observational data of plasma turbulence measured on December 9, 2016 (Phan et al. 2018; Stawarz et al. 2019). Indeed, one simulation setup that we include in this work is a 3D equivalent of the twodimensional (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 Sects. 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 Sect. 4.1 we give a precise definition of “current structure” in our simulations and define two nonlocal and overall shape factors called planarity 𝒫 and elongation ℰ to better highlight the 3D shape 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 𝒫 and elongation ℰ to characterize any such configuration. Then, in Sect. 4.2 we present three different methods by which it is possible to convert purely local measures into an estimate of ℰ and 𝒫 of current structures, and we show their effectiveness. In Sect. 5 we investigate the physical features of the current regions forming in the 3D simulations we make use of 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 9, 2016 (Phan et al. 2018; Stawarz et al. 2019). Finally, we summarize the results obtained and their importance in Sect. 6.
2. 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 the HVM model with kinetic ions and fluid neutralizing electrons 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 . As a result, the dimensionless electron skin depth is given by the electrontoion mass ratio, .
The ion distribution function f_{i} = f_{i}(x, v, t) evolves following the Vlasov equation that in dimensionless units reads as
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 electroninertia 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 lowfrequency regime). In Eq. (2) we assume an isothermal equation of state for the electron pressure, P_{e} = nT_{0e}, with a given initial electrontoion temperature ratio of T_{0e}/T_{0i} = 1. The initial ion distribution function is given by a Maxwellian distribution with a corresponding uniform temperature. Finally, the evolution of the magnetic field B is given by the following 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 a 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 (rms) amplitude of each component of the perturbed magnetic field, and the ratio between the rms of the magnetic perturbation and the equilibrium field. For all simulations, we sampled the velocity space using 51^{3} uniformly distributed grid points spanning [ − 5v_{th, i}, 5v_{th, i}] in each direction, where is the initial ion thermal velocity and β_{i} = 1. We set a reduced mass ratio of m_{i}/m_{e} = 100.
Simulation parameters.
In the following we refer to the 3D3V simulations (the first three listed in Table 1) using the following names: “SIM A”, “SIM Bwf” (weak forcing scenario), and “SIM Bsf” (strong forcing scenario). We also make use as a reference for comparison of a 2D3V hybrid VlasovMaxwell simulation, namely “SIM 2D” (last one in Table 1).
The reasons behind the choice of simulation parameters are discussed in the following section, also providing a brief overview of turbulence evolution and properties.
3. Context: Turbulent evolution and spectra
The simulations considered in this work are meant to reproduce the behavior of a β ≈ 1 turbulent plasma, typical of nearEarth’s environment and of the solar wind at about 1 AU (see for example, Chen 2016; Sahraoui et al. 2020 and references therein). As summarized in Table 1, each simulation initially injects fluctuations’ energy in a different wavenumber range, with different rms values of the initial magnetic fluctuations exploring a different number of subiongyroradius scales. In particular, SIM Bwf and SIM Bsf inject energy only slightly above the ion characteristic scales, and they are able to excite a wide range of kinetic scales before reaching their dissipation scale (that is, resolving about a decade of clean subion range turbulence before entering the dissipationdominated scales). In fact, as done in Califano et al. (2020) for the 2D3V case, the aim of both SIM Bwf and SIMsf is to mimic the conditions possibly encountered in the Earth’s magnetosheath, past the bow shock, where the occurrence of electrononly 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 2D 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 toward smaller and smaller scales. In the transition phase, (2) coherent structures, that is regions where the current density peaks (it usually happens at about one eddyturnover time) start to form. At the end, a fully developed turbulent state (3) is reached, with complex nonlinear dynamics with coherent structures continuously forming, merging, and/or being destroyed. We indicate the times at which the 3D simulations reach saturation (maximum of the rms of the current density) as t_{satA} (SIM A), t_{satB − wf} (SIM Bwf), and t_{satB − sf} (SIM Bsf); at these times, turbulence is fully developed (roughly speaking saturation times correspond to about 3 eddyturnover times). From now on, unless explicitly stated otherwise, the entire analysis will be carried out at saturation times.
In Fig. 1 we show the reduced onedimensional (1D), k_{∥}averaged, 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 and 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 Bwf and SIM Bsf 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 prevents from properly establishing critical balance (Schekochihin et al. 2009), this discrepancy is also observed in the spectrum emerging from the 3D case (SIM A, whose magneticfield spectrum seems to agree very well with the corresponding spectrum of SIM 2D). As also noticed in previous 3D3V simulations (Cerri et al. 2019), whether or not this feature in SIM A is due to the limited extent of the MHD range is unclear and it 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 Bsf; a hint of such a power law is also visible in SIM A, but the limited extent of the subion 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” kineticAlfvénwave cascade occurring at subion scales (Boldyrev & Perez 2012), and previously reported via 3D3V hybridVlasov simulations (Cerri et al. 2017a, 2018). The other two simulations, SIM 2D and SIM Bwf, exhibit a power law steeper than −8/3. A magneticfield spectrum close to in the kinetic range was already observed in previous 2D hybridkinetic simulations (see for example, Cerri et al. 2017b and references therein), sometimes attributed to the effect of a reconnectionmediated 2D cascade at subion scales (Cerri & Califano 2017; Franci et al. 2017), which, however, may continue to hold in three spatial dimensions as well (Loureiro & Boldyrev 2017; Mallet et al. 2017), although definitive evidence of such a 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 Bwf, on the other hand, a steeper spectrum may be attributed to a weaker cascade associated with the lower amplitude (with respect to SIM Bsf) of the injected fluctuations, rather than to dissipation effects associated with ionheating and/or fluctuations’ damping mechanisms in the β ≈ 1 regime considered here (see for example, Told et al. 2015; Sulem et al. 2016; Arzamasskiy et al. 2019 and references therein). In fact, in SIM Bwf, 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 ionheating mechanisms (and how they change between different simulations) is beyond the scope of this work and will be reported elsewhere.
Fig. 1. For our four simulations, 1D k_{∥}averaged magnetic energy spectra versus k_{⊥}, normalized so that they overlap at k_{⊥}d_{i} ≈ 2. 
In Fig. 2 we show for (a) SIM A and (b) SIM Bst, the 3D rendering of the current density magnitude in shaded isocontours; we would like to remind the reader that SIM Bst employs a simulation box size that, compared to the domain of SIM A, is about onethird in each direction. In blue, the regions of grid points which overcome a specific threshold for the current density J_{th}, defined in the next section, are shown.
Fig. 2. 3D rendering of regions with J ≥ J_{th} in blue for (a) SIM A and (b) SIM Bst. In gray, shadow color plots of the current density are shown. 
4. Current structures and magnetic configuration analysis
4.1. 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 such definitions, we also discuss the concept of “magnetic configuration” and the meaning of its shape, a concept which is thoroughly applied in the following sections.
Similarly to what was done in Uritsky et al. (2010) and Zhdankin et al. (2013), we adopt the formal definition of “current structure” as any connected region where the currentdensity magnitude, J ≡ J, is above a threshold value for all internal points, and the boundaries do not touch any other zone such as this. In particular, we stress that the specific definition of such a 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 , where and ⟨…⟩ denotes the spatial average over the whole simulation box. We want to remark that, by this procedure, any isolated point with an abovethreshold current density is not recognized as a current structure on its own, but it can be identified as part of some other current structure if the two are within a minimum distance of one gridpoint from each other. We call the “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 with ℓ_{max}). In this particular plane, we have a direction of strongest variation for the current density (associated with ℓ_{min}) and a direction of weakest variation (associated with ℓ_{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 being restricted 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 the impact is of slightly different definitions of ℓ_{min} on our analysis.
Once the “geometry” has been clarified, it is necessary to define the “shape” of the current structures as well. In a current structure, we call “shape” the quantity determined by two “structure shape factors” defined as follows:
These two parameters measure the tendency of the current structure to squeeze toward some elongated form (ℰ → 1) or to flatten (𝒫 → 1). In this way, we can represent each structure’s shape in the plane ℰ − 𝒫, and eventually classify shapes into categories. In particular, we adopt the following nomenclature for certain characteristic shapes: (i) “pseudo spheres” (low ℰ, low 𝒫), (ii) “cigars” (high ℰ, low 𝒫), (iii) “pancakes” (low ℰ, high 𝒫), and (iv) “knife blades” (high ℰ, high 𝒫). In Fig. 3 we show a schematic representation of a current structure with its characteristic lengths: ℓ_{min}, ℓ_{med}, and ℓ_{max}. We note that the three quantities of thickness, width, and length presented here are enough to define the shape of compact structures in 3D. In the case in which the holes would appear in current structures, new parameters should be added as a filling fraction or the hole dimensions. Such an upgrade will be investigated in future work.
Fig. 3. Schematic representation of a current structure with its characteristic lengths: ℓ_{min}, ℓ_{med}, and ℓ_{max}. 
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}, and σ_{min}. The quantities , , and 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, and not directly comparable, due to the presence of normalization factors. What is really significant is the ratio between eigenvalues, as we see in the next section. Indeed, by 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 , , and in the exact same way we did with characteristic lengths and to 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 owing to a current structure. On the contrary, the “structure shape factors” provide an overall and nonlocal picture of a current structure. As ℰ and 𝒫 do for current structures, ℰ and 𝒫 do measure the tendency of magnetic configurations toward elongation (if ℰ → 1) or squashing into sheets (if 𝒫 → 1). Similarly, the position of a configuration in the ℰ − 𝒫 plane allows one 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 only focused on satellite data and thus on local measurements, providing a pointbypoint value for ℰ and 𝒫 along spacecrafts’ trajectories. On the contrary here, taking advantage of the 3D data of a simulation, we can also give an overall and nonlocal picture of current structures (via ℰ and 𝒫).
4.2. The shape of current structures
The first goal of this paper is to determine whether it is possible to estimate the “nonlocal” overall shape of a current structure by local measurements, that is 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 ℰ and 𝒫 as given by our “nonlocal” overall method in Eqs. (4) and (5) to be used as “reference values” for the estimates obtained via the following local methods. Firstly, “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 follows:
where the subscript “cen” indicates that the values are calculated at the structure’s center. We note that this method closely resembles the one adopted by Servidio et al. (2009). Secondly, “NB” method. This approach aims at inferring the “structure shape factors”, ℰ and 𝒫, from the “magneticconfiguration shape factors”, ℰ and 𝒫, 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}, and σ_{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 is representative for the whole structure. However, to easily make a comparison with data for satellites and/or to increase the dataset analyzed, one can apply this method even to points different from the central one. Since the magnetic field and current density are strongly related (in our model, the current density is the curl of the magnetic field), we also expect ℰ^{NB} and 𝒫^{NB} to resemble ℰ^{HJ} and 𝒫^{HJ}. Finally, “AV” method. This approach considers an average of the magnetic configurations occurring inside each structures, that is to say shape parameters are estimated as follows:
where ⟨…⟩_{str} is the average over all points belonging to a single current structure. We note that this procedure resembles the NB method, but it can be carried out without knowing where the center of a current structure is located.
4.3. Method comparison
Our goal here is to verify the accuracy of the different methods HJ, NB, and AV at estimating ℰ and 𝒫, which are obtained by looking at the overall and nonlocal shape of the current structure and are thus considered our reference values. To establish how accurate the different methods detailed in Sect. 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 Fig. 4 we show how ℰ and 𝒫 are distributed when defined through ℓ_{max}, ℓ_{med}, and ℓ_{min} (Eqs. (4) and (5)), while in Fig. 5 we report the results of the estimates from the HJ, NB, and AV methods (that is, distributions for ℰ^{HJ} and 𝒫^{HJ}, ℰ^{NB} and 𝒫^{NB}, and ℰ^{AV} and 𝒫^{AV}, respectively, alongside their ratios to the reference “nonlocal” overall ℰ and 𝒫 values for each structure).
Fig. 4. Occurrence distributions of elongation ℰ versus planarity 𝒫 (as computed using the “nonlocal” overall method discussed in Sect. 4.1) for SIM A. 
The physical picture that emerges from Fig. 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 the mean, median, and standard deviation, see the first two rows of Table 2), suggesting “knifeblade” 3D configurations. The results from HJ and AV methods are in very good agreement with the overall occurrence distribution of ℰ and 𝒫, as shown in the top row panels of Fig. 5. For what concerns method NB, it gives the same physical picture of the other local and nonlocal 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 toward the upperright 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 Fig. 5. More specifically, Fig. 5d shows that the distribution of ℰ/ℰ^{HJ} is peaked around 0.9, and that 𝒫/𝒫^{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 the HJ method (for the mean, median, and standard deviation see the third and fourth rows in Table 2). The same agreement is shown in Fig. 5e between the control values and the NB method estimates, as shown by the mean, median, and standard deviation reported in Table 2 (fifth and sixth rows). Finally, Fig. 5f 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 the seventh and eighth row of Table 2.
Fig. 5. Occurrence distributions of elongation versus planarity, for SIM A. Top row panels: (a) HJ method, (b) NB method, and (c) AV method. Bottom row panels: (d) ℰ/ℰ^{HJ} versus 𝒫/𝒫^{HJ}, (e) ℰ/ℰ^{NB} versus 𝒫/𝒫^{NB}, and (f) ℰ/ℰ^{AV} versus 𝒫/𝒫^{AV}, thus the comparison between our local methods and our overall nonlocal reference values ℰ and 𝒫. The x and y ranges in panels d–f have been adapted in order to produce a better visualization. In doing so, some cases have been excluded including the following: 21 out of the total 1043 current structures (CSs) examined in panel d, 23 out of 1043 in panel e, and two out of 1043 in panel f. 
Mean, median, and standard deviation for method comparison.
Despite the overall picture providing a very good agreement between reference and local methods, we note the presence of a slightly, yet 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 a proper current threshold is defined.
5. 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 investigate the magnetic configuration shape factors (planarity and elongation), the characteristic eigenvalues (proportional to the characteristic scale lengths), and the orientation of the magnetic field and current density. We use the N tensor and its eigenvalues, but without restricting our analysis to points owing to the current structures, as for NB and AV methods, thus we essentially apply the MCA method as first described 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 consider separately (1) all points belonging to current structures and (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 a satellite along 1D trajectories. Indeed, in the uniform sampling we picked one point for every three in each direction. Because the typical dimensions of a current structure are, on average, big enough to include more than three grid points in at least two directions, for each structure several points are collected, as is the case for a continuous satellite sampling. Investigating the possible difference between a 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 3D3V simulations with an analogous analysis applied to two intervals of highresolution (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 threestep selection on the abovementioned MMS data has been applied so that the magnetic configuration analysis has been performed on only 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 the original data are kept; for details, see Appendix B).
5.1. Analysis of shape factors for magnetic configurations
In Fig. 6 we show the occurrence distribution of ℰ versus 𝒫 in our set of 3D simulations, namely SIM A, SIM Bwf, and SIM Bsf (left, center, and right column, respectively). This can be seen 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) and in the bottom row only for points belonging to what is identified as a “current structure” (that is, those regions where the current density is above J_{th}).
Fig. 6. Occurrence distribution of elongation ℰ versus planarity 𝒫 for the three different simulations (three columns), in the top row for generic points and in the bottom row for the points which belong to current structures. 
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 Fig. 6 with those obtained in Sect. 4.3 for the same simulation (that is, see Figs. 4 and 5a–c). We further point out that for SIM Bwf, this discrepancy between the two distributions (viz., generic points versus only points belonging to current overdensities) seems to be particularly emphasized. In the following we 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 Fig. 7 where we show the 1D normalized occurrence distribution of ℰ and 𝒫 for the three simulations SIM A, SIM Bwf, and SIM Bsf for generic points (magenta line) and for points belonging to current structures (black line). We have superposed these distributions obtained from simulations to the ones obtained using MMS satellite data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019) (see Appendix B for details).
Fig. 7. Normalized occurrence distributions of 𝒫 and ℰ for the three simulations SIM A, SIM Bwf, and SIM Bsf, distinguishing between statistics on generic points (magenta line) and points which belong to current structures (black line). We superpose these distributions, obtained from simulations, to the ones obtained using satellite data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. 
As in the colormaps of Fig. 6, the histograms in Fig. 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 𝒫 distributions less skewed toward 𝒫 ≈ 1, with a peak around 0.6 − 0.8 (it is important to note the log scale for the y axis). On the contrary, the normalized histogram of 𝒫 for points owing to current structures peaks around 0.8 − 1.0. For all three simulations, the values for the mean, median, and standard deviation are reported in Table 3 (first row for generic points and second 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 the mean, median, and standard deviation in Table 3 (third row for generic points and fourth row for current structures). In summary, the emerging picture indicates a majority of “bladelike” 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 a high planarity confirming the intuitive picture of nearly 2D current sheets. On the contrary, for generic points which are picked up from almost everywhere, including from regions where the current density is low, the magnetic configuration has a different behavior.
Mean, median, and standard deviation for configuration analysis.
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 Bsf. 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 subset of “generic points” from our simulations. Second, simulation SIM Bsf has a setup which is the 3D equivalent of the one employed in Califano et al. (2020), which has already been proven 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 Bsf and these particular MMS intervals.
5.2. Analysis of eigenvalues’ occurrence distribution and characteristic scale lengths in simulations and satellite data
In Fig. 8 we show the normalized occurrence distribution of N’s eigenvalues: (a) σ_{min}, (b) σ_{med}, and (c) σ_{max} for SIM A, SIM Bwf, and SIM Bsf, respectively. We recall that these eigenvalues are interpreted as proportional to the squared inverse of the local variation length (see Sect. 4.1). We superpose these distributions, obtained from simulations, to the ones obtained using satellites data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019) (with the additional selection discussed at the beginning of this section; see Appendix B for details). In the following, we note that the MMS eigenvalues are normalized using the local value of d_{i}.
Fig. 8. Normalized occurrence distribution of N’s eigenvalues: (a) σ_{min}, (b) σ_{med}, and (c) σ_{max} for SIM A, SIM Bwf, and SIM Bsf, respectively. In magenta distributions are shown for generic points, and in black they are shown for points belonging to current structures. We superposed these distributions, obtained from simulations, to the ones obtained using satellites data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. We note the logarithmic scale on the x axis. 
Let us first consider the occurrence distribution of the smallest eigenvalue, which is related to the largest characteristic length of the magnetic configuration. This is shown in Col. a of Fig. 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 significantly differ from that of a structure where the current peaks. Now let us instead consider the occurrence distributions for the other two eigenvalues (related to the median and shortest variation length of the local magnetic configurations, shown in Cols. b and c of Fig. 8, respectively). In all the simulations, the occurrence distributions of these eigenvalues 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, that is the one associated with 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 toward 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 Bwf and SIM Bsf, 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 Bsf, rather than for SIM Bwf. Thus, simulation SIM Bsf develops current structures with a typical scale length smaller than those of the other simulations. This is also true for generic points, even if less pronounced. This type of feature is likely a consequence of the fact that the two simulations develop different turbulent regimes (see discussion in Sect. 3). In fact, by injecting largeramplitude magnetic fluctuations in SIM Bsf than in SIM Bwf, a shallower subionscale spectrum develops and, consequently, a larger amount of turbulent power reaches the smallest scales (cf. the magneticfield spectrum in Fig. 1). Whether this is actually due to a weaklike versus stronglike turbulent regime (that is, if it could be a consequence of whether critical balance and/or dynamic alignment establishes or not) or to other effects (for example, development of a significant amount of electrononly 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, and standard deviation values in the fifth and eighth rows in Table 3, comparing SIM Bwf and SIM Bsf quantitatively.
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 Sect. 5.1, this very good agreement with SIM Bsf was somewhat expected since a similar simulation setup in the 2D case has already been proven capable of qualitatively reproducing the turbulent and reconnection regime observed in these intervals (Califano et al. 2020). Thus, we proved once again that a relatively largeamplitude injection close to the ionkinetic 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 a setup constitutes a quite interesting hint, the actual mechanism that could be behind such a peculiar injection (for example, the bow shock itself or the subsequent occurrence of microinstabilities) is still unclear and requires further investigations, especially from the observational point of view.
5.3. Analysis of the orientation of the 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 by restricting analysis to only points which belong to current structures.
In Fig. 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 these distributions to those obtained by using satellites data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019) (see Appendix B for details).
Fig. 9. Normalized occurrence distribution of b ⋅ e_{min} (left column) and j ⋅ e_{min} (right column) for all three simulations. Magenta is used to show the distribution for generic points, and black corresponds to points belonging to current structures. The distribution refers to times at fully developed turbulence, in particular t_{satA} for SIM A, t_{satB − wf} for SIM Bwf, and t_{satB − sf} for SIM Bsf. The distribution refers to times at fully developed turbulence, in particular t_{satA} for SIM A, t_{satB − wf} for SIM Bwf, and t_{satB − sf} for SIM Bsf. We superposed these distributions, obtained from simulations, to the ones obtained using satellite data (dotted blue line) from the two highresolution 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 the 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 the current and magnetic field, which is typical of smallscale turbulent structures (see for example, 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 Bwf rather than with SIM Bsf. 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 Bwf, 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 Figs. 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.
5.4. Generic points or “background”
We briefly explain what it 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 a simulation (which is instead the subset 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 subset of “generic” points has been created by considering one point for every 27 (corresponding to a collection of ∼10^{6} points), that is around ∼4% of the original set. Analogously, the total number of points which overcome the threshold on current density at time t_{satA} are ∼1.1 × 10^{6} (that is, roughly 4% of the total). Moreover, in our subset of “generic” points those which overcome the threshold on current density are 4.8 × 10^{4}, again roughly 4% of the subset. Thus, the subset of “generic” points that is being considered should adequately represent the ensemble of points of our simulation box. As recently discussed, the points which constitute the current structures are always a small percentage of the total number of points in a set or subset (that is, the filling factor of the current structures in a volume is usually very small). Therefore, we might wonder if the behavior of such a small subgroup can 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 Fig. 10, which shows the occurrence distribution of elongation ℰ versus planarity 𝒫 for (a) only those “generic” points that do not belong to current structures (roughly 4%) and (b) the whole subset of “generic” points (that is, including those belonging to current structures). In fact, there is no noticeable difference between the two distributions. In particular, for both cases, the elongation and planarity have the same mean, median, and standard deviation (see Table 3, rows 9 to 12).
Fig. 10. For t_{satA} SIM A, the occurrence distribution of elongation versus planarity for (a) “generic” points excluding the ones which belong to current structures (roughly 4%) and (b) “generic” points. 
The behavior of the current structures cannot emerge if we consider “generic” points, which are the equivalent of a long continuous sample in satellite data, without introducing any additional selection on the turbulent time series. Therefore, in order to systematically study the statistical behavior of the shape of current structures in satellites data, we suggest fixing a threshold on the current density and only performing analysis on those regions that overcome this threshold. This procedure is similar to what has been done when applying the AV method discussed in Sect. 4.2 to our numerical simulations. A similar kind of suggestion (to use a threshold to identify strong current points) has also been proposed in other works such as Bruno et al. (2001), Greco et al. (2009), Osman et al. (2012) and SorrisoValvo et al. (2018).
6. Conclusions
In this work we have conducted a widespanning analysis of overall nonlocal current structures and local magnetic configurations emerging in plasma turbulence. The analysis is based on three 3D3V HVM simulations with different box sizes, resolution, and energy injection scales. We first focused on the characterization of the shape of current structures. Our “reference method” defines two parameters, elongation ℰ and planarity 𝒫, 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 by also employing three different local methods, namely HJ, NB, and AV (see Sect. 4.2 for details on the methods). The rigorous computation of ℰ and 𝒫 via the “reference” nonlocal method can only be performed on simulation data. The HJ and NB methods can also only be performed on simulation data since they require being 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 (that is, the “nonlocal” and the AV methods), since they require being 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 also be applied with minor modifications 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 (that is, “reference and nonlocal”, HJ, NB, and AV), we have analyzed the distribution of shape factors (that is, planarity and elongation) for the emerging current structures, with the result that all methods coherently find that they are composed of mainly “knifeblade”like structures. This picture is different from the one that emerges in Meyrand & Galtier (2013) where they claim that there is a presence of mostly cigarlike structures (they use the expression “filamentlike”) through visual inspection of their 3D HallMHD simulations of turbulence. This suggests that ionkinetic effects (that is, beyond just the Hall term) and/or electroninertia terms could significantly affect (and likely be required to correctly describe) the development of a current structure in plasma turbulence across the transition range and at subion 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 magneticfield and currentdensity 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 (that is, including, but not limited to, structures). In particular, the main statistical difference between generic points and those located inside a current structure is in the results obtained for the distribution of planarity 𝒫: “knifeblade” 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 𝒫 and elongation ℰ is coherent with the distribution we obtain if we apply the same kind of analysis to highresolution (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 Bwf and SIM Bsf, for which all the three characteristic length scales are shorter in the strongforcing scenario (SIM Bsf) than they are when loweramplitude fluctuations are injected (SIM Bwf). We remind the reader that in both simulations SIM Bwf and SIM Bsf, the energy is injected only slightly above the ion characteristic scales, as done in Califano et al. (2020), which was able to reproduce, in a simplified 2D3V configuration, the electrononly 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 Bsf, thus suggesting that similar, smallscale current structures require a high level of forcing acting close to the ionkinetic 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 the MMS of Phan et al. (2018) and Stawarz et al. (2019) are very good. In particular, we note a good agreement with the distributions for generic points in SIM Bwf.
The entire analysis presented in this work clearly highlights how magnetic configurations inside a current structure exhibit peculiar features that can only be retrieved by solely considering those points where J attains values above a certain threshold. Indeed, we have shown in Sect. 5.4 that the behavior of the current structures cannot emerge if we consider “generic” points, which is 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 is only a small percentage of the total. Therefore, in order to systematically study the shape of current structures in satellites data, we suggest fixing a proper threshold on the current density and consequently considering only regions that overcome this threshold for the subsequent analysis. Even better, one could apply the AV method that we have described in Sect. 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 not only be useful for the analysis of turbulence simulations, but also for observative studies. Indeed, (1) we have validated the possibility of applying local methods, which are the only ones applicable on satellite data, to infer the overall nonlocal shape of current structures; (2) we conjecture that imposing a proper threshold on the current density would be beneficial 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 exclusively considering points within current structures with respect to what emerges from the points belonging to the rest of the turbulent environment; and (4) we have shown via such a magnetic configuration analysis that, when there is a mechanism that injects relatively highlevel fluctuations close to the ionkinetic scales, our 3D3V simulation can reproduce the structures that emerge in MMS data for the periods studied by Phan et al. (2018) and Stawarz et al. (2019).
Acknowledgments
This work has received funding from the European Unions Horizon 2020 research and innovation programme under grant agreement No. 776262 (AIDA, www.aidaspace.eu). We acknowledge the CINECA ISCRA initiative and the EU PRACE initiative (grant n. 2017174107) for awarding us access to the supercomputer Marconi, CINECA, Italy, where the calculations were performed. FC thanks Dr. M. Guarrasi (CINECA, Italy) for his essential contribution for his contribution on code implementation on Marconi. S.S.C. was supported by NASA Grant No. NNX16AK09G issued through the Heliophysics Supporting Research Program, and by the MaxPlanck/Princeton Center for Plasma Physics (NSF grant PHY1804048).
References
 Agudelo Rueda, J. A., Verscharen, D., Wicks, R. T., et al. 2021, J. Plasma Phys., 87, 905870228 [CrossRef] [Google Scholar]
 Arzamasskiy, L., Kunz, M. W., Chandran, B. D. G., & Quataert, E. 2019, ApJ, 879, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Boldyrev, S., & Perez, J. C. 2012, ApJ, 758, L44 [NASA ADS] [CrossRef] [Google Scholar]
 Bruno, R., Carbone, V., Veltri, P., Pietropaolo, E., & Bavassano, B. 2001, Planet. Space Sci., 49, 1201 [Google Scholar]
 Burch, J. L., Moore, T. E., Torbert, R. B., & Giles, B. L. 2016, Space Sci. Rev., 199, 5 [CrossRef] [Google Scholar]
 Califano, F., Cerri, S. S., Faganello, M., et al. 2020, Front. Phys., 8, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Camporeale, E., SorrisoValvo, L., Califano, F., & Retinò, A. 2018, Phys. Rev. Lett., 120, 125101 [CrossRef] [Google Scholar]
 Carbone, V., Veltri, P., & Mangeney, A. 1990, Phys. Fluids A: Fluid Dyn., 2, 1487 [NASA ADS] [CrossRef] [Google Scholar]
 Cerri, S. S., & Califano, F. 2017, New J. Phys., 19, 025007 [NASA ADS] [CrossRef] [Google Scholar]
 Cerri, S. S., Servidio, S., & Califano, F. 2017a, ApJ, 846, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Cerri, S. S., Franci, L., Califano, F., Landi, S., & Hellinger, P. 2017b, J. Plasma Phys., 83, 705830202 [CrossRef] [Google Scholar]
 Cerri, S. S., Kunz, M. W., & Califano, F. 2018, ApJ, 856, L13 [CrossRef] [Google Scholar]
 Cerri, S. S., Grošelj, D., & Franci, L. 2019, Front. Astron. Space Sci., 6, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C. H. K. 2016, J. Plasma Phys., 82, 535820602 [Google Scholar]
 Comisso, L., & Sironi, L. 2019, ApJ, 886, 122 [Google Scholar]
 De Giorgio, E., Servidio, S., & Veltri, P. 2017, Sci. Rep., 7, 13849 [NASA ADS] [CrossRef] [Google Scholar]
 Dong, C., Wang, L., Huang, Y.M., Comisso, L., & Bhattacharjee, A. 2018, Phys. Rev. Lett., 121, 165101 [NASA ADS] [CrossRef] [Google Scholar]
 Dunlop, M., Southwood, D., Glassmeier, K.H., & Neubauer, F. 1988, Adv. Space Res., 8, 273 [CrossRef] [Google Scholar]
 Fadanelli, S., Lavraud, B., Califano, F., et al. 2019, J. Geophys. Res.: Space Phys., 124, 6850 [NASA ADS] [CrossRef] [Google Scholar]
 Franci, L., Cerri, S. S., Califano, F., et al. 2017, ApJ, 850, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Gingell, I., Schwartz, S. J., Eastwood, J. P., et al. 2020, J. Geophys. Res.: Space Phys., 125, e27119 [NASA ADS] [CrossRef] [Google Scholar]
 Gosling, J. T., & Phan, T. D. 2013, ApJ, 763, L39 [Google Scholar]
 Greco, A., Matthaeus, W. H., Servidio, S., & Dmitruk, P. 2009, Phys. Rev. E, 80, 046401 [NASA ADS] [CrossRef] [Google Scholar]
 Greco, A., Matthaeus, W. H., Perri, S., et al. 2018, Space Sci. Rev., 214, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Grošelj, D., Cerri, S. S., Bañón Navarro, A., et al. 2017, ApJ, 847, 28 [CrossRef] [Google Scholar]
 Khabarova, O., Malandraki, O., Malova, H., et al. 2021, Space Sci. Rev., 217, 38 [CrossRef] [Google Scholar]
 Loureiro, N. F., & Boldyrev, S. 2017, ApJ, 850, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2017, J. Plasma Phys., 83, 905830609 [NASA ADS] [CrossRef] [Google Scholar]
 Matthaeus, W. H., Yang, Y., Wan, M., et al. 2020, ApJ, 891, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Meyrand, R., & Galtier, S. 2013, Phys. Rev. Lett., 111, 264501 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, A. B., Teaca, B., Told, D., et al. 2016, Phys. Rev. Lett., 117, 245101 [NASA ADS] [CrossRef] [Google Scholar]
 Osman, K. T., Matthaeus, W. H., Hnat, B., & Chapman, S. C. 2012, Phys. Rev. Lett., 108, 261103 [NASA ADS] [CrossRef] [Google Scholar]
 Osman, K. T., Matthaeus, W. H., Gosling, J. T., et al. 2014, Phys. Rev. Lett., 112, 215002 [NASA ADS] [CrossRef] [Google Scholar]
 Papini, E., Franci, L., Landi, S., et al. 2019, ApJ, 870, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Passot, T., Henri, P., Laveder, D., & Sulem, P.L. 2014, Eur. Phys. J. D, 68, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Pecora, F., Greco, A., Hu, Q., et al. 2019, ApJ, 881, L11 [Google Scholar]
 Phan, T. D., Eastwood, J. P., Shay, M. A., et al. 2018, Nature, 557, 202 [Google Scholar]
 Podesta, J. J. 2017, J. Geophys. Res.: Space Phys., 122, 2795 [NASA ADS] [CrossRef] [Google Scholar]
 Pollock, C., Moore, T., Jacques, A., et al. 2016, Space Sci. Rev., 199, 331 [CrossRef] [Google Scholar]
 Russell, A. J. B., Yeates, A. R., & Eastwood, J. P. 2015, Astron. Geophys., 56, 6.18 [NASA ADS] [CrossRef] [Google Scholar]
 Sahraoui, F., Hadid, L., & Huang, S. 2020, Rev. Mod. Plasma Phys., 4, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Cowley, S. C., Dorland, W., et al. 2009, ApJS, 182, 310 [NASA ADS] [CrossRef] [Google Scholar]
 Servidio, S., Matthaeus, W. H., Shay, M. A., Cassak, P. A., & Dmitruk, P. 2009, Phys. Rev. Lett., 102, 115003 [Google Scholar]
 Shen, C., Dunlop, M., Li, X., et al. 2007, J. Geophys. Res.: Space Phys., 112 [Google Scholar]
 Shi, Q. Q., Shen, C., Pu, Z. Y., et al. 2005, Geophys. Res. Lett., 32 [Google Scholar]
 Shi, Q. Q., Tian, A. M., Bai, S. C., et al. 2019, Space Sci. Rev., 215, 35 [CrossRef] [Google Scholar]
 SorrisoValvo, L., Carbone, F., Perri, S., et al. 2018, Sol. Phys., 293, A10 [Google Scholar]
 Stawarz, J. E., Eastwood, J. P., Phan, T. D., et al. 2019, ApJ, 877, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Sulem, P. L., Passot, T., Laveder, D., & Borgogno, D. 2016, ApJ, 818, 66 [CrossRef] [Google Scholar]
 TenBarge, J. M., & Howes, G. G. 2013, ApJ, 771, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Told, D., Jenko, F., TenBarge, J. M., Howes, G. G., & Hammett, G. W. 2015, Phys. Rev. Lett., 115, 025003 [NASA ADS] [CrossRef] [Google Scholar]
 Uritsky, V., Pouquet, A., Rosenberg, D., Mininni, P., & Donovan, E. 2010, Phys. Rev. E, 82, 056326 [NASA ADS] [CrossRef] [Google Scholar]
 Valentini, F., Trávníček, P., Califano, F., Hellinger, P., & Mangeney, A. 2007, J. Comput. Phys., 225, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Vech, D., Mallet, A., Klein, K. G., & Kasper, J. C. 2018, ApJ, 855, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Zhdankin, V., Uzdensky, D. A., Perez, J. C., & Boldyrev, S. 2013, ApJ, 771, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Zhdankin, V., Boldyrev, S., Perez, J. C., & Tobias, S. M. 2014, ApJ, 795, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Zhdankin, V., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2017, Phys. Rev. Lett., 118, 055103 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Alternative definitions for structure’s dimensions
In our work, we define 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. Here, we present 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 − turnover2D}: in red the estimations given by the method used in the main text are shown (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 are shown (full width at half maximum of the current density profile, called “FWHM method”). We chose the 2D simulations for this comparison and the eddyturnover time rather than the saturation time just because in this simulation, at this time, there are not many current structures, thus it is easier to interpret the scatter plot. The current density was 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 for the current density. In the subpanel, 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}, the thickness estimated using the “FWHM” method is generally bigger than the one estimated using the “mask method”. On the contrary, when the current density is less 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 ℰ and 𝒫 which emerges, which is a bit different from the one of Figure 4. In particular, to quantify, now we have the following: the elongation mean ∼0.8, the median ∼0.9, and the standard deviation ∼0.2, and instead the planarity mean ∼0.4, the median ∼0.4, and the standard deviation ∼0.2. The distribution is in less agreement with the ones for the HJ, NB, and AV methods shown in Figure 5, and most of the structures have a planarity in between 0.20.6, thus they are more filamentlike.
Fig. A.1. Scatter plot thickness versus current density for SIM 2D at t_{eddy − turnover2D} to show the differences between the two methods in estimating the thickness. The current density was computed at the center of each current structure. By “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, by “FWHM method”, we refer to the alternative method explained in Appendix A which computes the thickness using full width at half maximum of the interpolated profile of J. In the subpanel we show a schematic representation to explain the difference in estimations between the “FWHM method” and the “mask method”. 
Fig. A.2. Occurrence distributions of elongation ℰ versus planarity 𝒫 at t = t_{satA} for SIM A. Planarity here was computed using the thickness given by the alternative method explained into Appendix A. 
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 burstresolution 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. The choice of these MMS data intervals is motivated by the fact that our simulation setup (in particular for SIM Bwt and SIM Bst) 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 have been taken from the FluxGate Magnetometer (FGM  see Russell et al. 2015), while the Fast Plasma Investigation (FPI  Pollock et al. 2016) has provided the density and pressure measures which contributed to determining 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 fleets. 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), that is by performing a linear estimation of the magnetic field gradient, then combining the resulting values with the fourspacecraft averaged magnetic field data.
In order to compare results from MCA applied over simulation and MMS data, we needed to apply three different levels of selection to the latter. In particular, (1) we selected points for which at least two of the eigenvalues of the N matrix are welldetermined by MMS measurements; (2) we selected data with β ∈ [0.3, 3]; and (3) we selected data with a resolution comparable with the one of our simulations.
Furthermore, the first selection was performed exactly as in Fadanelli et al. (2019), that is by calculating the average interspacecraft distance ℓ_{SC} at each instant and then at setting a minimal resolution threshold for the eigenvalues of N, being understood that any eigenvalue below such a threshold is not wellresolved. Requiring that at least two of the eigenvalues are wellresolved signifies that, on the one hand, we accept uncertainty on the least eigenvalue and elongation measures but, on the other hand, we are still able to determine the direction of minimum variation correctly (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 the smallest eigenvalue correctly, which is generally extremely small and therefore easily falls below the MMS resolution threshold. Choosing the two welldetermined eigenvalue selection criterion for the two turbulence intervals that we consider here implies that over 99% of original data are retained by this first procedure.
The second selection just described was performed over the ion plasma beta obtained by the ratio of spacecraftaveraged 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 available 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 the precision of MMS data and of simulations, we introduce two quantities that we call “resolution factors” and they are defined as follows. For the MMS (satellite) measurements, the resolution factor is the interspacecraft separation divided by the ion inertial length (we note that the only possible derivative calculation in this case follows from a linear interpolation of satellite measures). For the simulations, on the other hand, the resolution factor is 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 Bwf and SIM Bsf, 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 200000 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.
All Tables
All Figures
Fig. 1. For our four simulations, 1D k_{∥}averaged magnetic energy spectra versus k_{⊥}, normalized so that they overlap at k_{⊥}d_{i} ≈ 2. 

In the text 
Fig. 2. 3D rendering of regions with J ≥ J_{th} in blue for (a) SIM A and (b) SIM Bst. In gray, shadow color plots of the current density are shown. 

In the text 
Fig. 3. Schematic representation of a current structure with its characteristic lengths: ℓ_{min}, ℓ_{med}, and ℓ_{max}. 

In the text 
Fig. 4. Occurrence distributions of elongation ℰ versus planarity 𝒫 (as computed using the “nonlocal” overall method discussed in Sect. 4.1) for SIM A. 

In the text 
Fig. 5. Occurrence distributions of elongation versus planarity, for SIM A. Top row panels: (a) HJ method, (b) NB method, and (c) AV method. Bottom row panels: (d) ℰ/ℰ^{HJ} versus 𝒫/𝒫^{HJ}, (e) ℰ/ℰ^{NB} versus 𝒫/𝒫^{NB}, and (f) ℰ/ℰ^{AV} versus 𝒫/𝒫^{AV}, thus the comparison between our local methods and our overall nonlocal reference values ℰ and 𝒫. The x and y ranges in panels d–f have been adapted in order to produce a better visualization. In doing so, some cases have been excluded including the following: 21 out of the total 1043 current structures (CSs) examined in panel d, 23 out of 1043 in panel e, and two out of 1043 in panel f. 

In the text 
Fig. 6. Occurrence distribution of elongation ℰ versus planarity 𝒫 for the three different simulations (three columns), in the top row for generic points and in the bottom row for the points which belong to current structures. 

In the text 
Fig. 7. Normalized occurrence distributions of 𝒫 and ℰ for the three simulations SIM A, SIM Bwf, and SIM Bsf, distinguishing between statistics on generic points (magenta line) and points which belong to current structures (black line). We superpose these distributions, obtained from simulations, to the ones obtained using satellite data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. 

In the text 
Fig. 8. Normalized occurrence distribution of N’s eigenvalues: (a) σ_{min}, (b) σ_{med}, and (c) σ_{max} for SIM A, SIM Bwf, and SIM Bsf, respectively. In magenta distributions are shown for generic points, and in black they are shown for points belonging to current structures. We superposed these distributions, obtained from simulations, to the ones obtained using satellites data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. We note the logarithmic scale on the x axis. 

In the text 
Fig. 9. Normalized occurrence distribution of b ⋅ e_{min} (left column) and j ⋅ e_{min} (right column) for all three simulations. Magenta is used to show the distribution for generic points, and black corresponds to points belonging to current structures. The distribution refers to times at fully developed turbulence, in particular t_{satA} for SIM A, t_{satB − wf} for SIM Bwf, and t_{satB − sf} for SIM Bsf. The distribution refers to times at fully developed turbulence, in particular t_{satA} for SIM A, t_{satB − wf} for SIM Bwf, and t_{satB − sf} for SIM Bsf. We superposed these distributions, obtained from simulations, to the ones obtained using satellite data (dotted blue line) from the two highresolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details. 

In the text 
Fig. 10. For t_{satA} SIM A, the occurrence distribution of elongation versus planarity for (a) “generic” points excluding the ones which belong to current structures (roughly 4%) and (b) “generic” points. 

In the text 
Fig. A.1. Scatter plot thickness versus current density for SIM 2D at t_{eddy − turnover2D} to show the differences between the two methods in estimating the thickness. The current density was computed at the center of each current structure. By “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, by “FWHM method”, we refer to the alternative method explained in Appendix A which computes the thickness using full width at half maximum of the interpolated profile of J. In the subpanel we show a schematic representation to explain the difference in estimations between the “FWHM method” and the “mask method”. 

In the text 
Fig. A.2. Occurrence distributions of elongation ℰ versus planarity 𝒫 at t = t_{satA} for SIM A. Planarity here was computed using the thickness given by the alternative method explained into Appendix A. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.