© M. Sisti et al. 2021

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 near-Earth 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 scale-to-scale 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) 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 a definition of characteristic lengths was employed, for instance, by Zhdankin et al. (2013) to characterize the current structures emerging in “reduced-MHD” (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 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 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 multi-spacecraft 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 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 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 three-dimensional 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 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 two-dimensional (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 $d_{\mathrm{i}}={\mathit{v}}_{\mathrm{A}}{\mathrm{\Omega }}_{\mathrm{ci}}^{-1}$. As a result, the dimensionless electron skin depth is given by the electron-to-ion mass ratio, $d_{\mathrm{e}}=\sqrt{m_{\mathrm{e}}/m_{\mathrm{i}}}$.

The ion distribution function f_i = f_i(x, v, t) evolves following the Vlasov equation that in dimensionless units reads as

$$\begin{array}{c}\hfill \frac{\partial f_{\mathrm{i}}}{\partial t}+\mathit{v}·\mathrm{\nabla }{f}_{\mathrm{i}}+(\mathit{E}+\mathit{v}\times \mathit{B})·\frac{\partial f_{\mathrm{i}}}{\partial \mathit{v}}=0,\end{array}$$(1)

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:

$$\begin{array}{c}\hfill (1-d_{\mathrm{e}}^2{\mathrm{\nabla }}^2)\mathit{E}=-{\mathit{u}}_{\mathrm{i}}\times \mathit{B}+\frac{\mathit{J}\times \mathit{B}}{n}-\frac{\mathrm{\nabla }{P}_{\mathrm{e}}}{n}+\frac{d_{\mathrm{e}}^2}{n}\mathrm{\nabla }·({\mathit{u}}_{\mathrm{i}}\mathit{J}+\mathit{J}{\mathit{u}}_{\mathrm{i}}-\frac{\mathit{JJ}}{n})·\end{array}$$(2)

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 Eq. (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 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:

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=-\mathrm{\nabla }\times \mathit{E}.\end{array}$$(3)

We take an initial magnetic field given by a uniform background field, B₀e_z where B₀ = 1, with a superimposed small-amplitude 3D perturbation, δB = δB_xe_x + δB_ye_y + δB_ze_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 wave-vector 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³ uniformly distributed grid points spanning [ − 5v_th, i, 5v_th, i] in each direction, where ${\mathit{v}}_{\mathrm{th},\mathrm{i}}=\sqrt{{\beta }_i/2}{\mathit{v}}_{\mathrm{A}}$ is the initial ion thermal velocity and β_i = 1. We set a reduced mass ratio of m_i/m_e = 100.

In the following we 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, 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 near-Earth’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 sub-ion-gyroradius scales. In particular, SIM B-wf and SIM B-sf 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 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 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 eddy-turnover 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 B-wf), and t_{satB − sf} (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, the entire analysis will be carried out at saturation times.

In Fig. 1 we show the reduced one-dimensional (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 $k_{\perp }^{-5/3}$ and $k_{\perp }^{-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 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 magnetic-field spectrum seems to agree very well with the corresponding spectrum of SIM 2D). As also noticed in previous 3D-3V 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 B-sf; a hint of such a power law is also visible 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 & Perez 2012), and previously reported via 3D-3V hybrid-Vlasov simulations (Cerri et al. 2017a, 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_{\perp }^{-3}$ in the kinetic range was already observed in previous 2D hybrid-kinetic simulations (see for example, Cerri et al. 2017b and references therein), sometimes attributed to the effect of a reconnection-mediated 2D cascade at sub-ion 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 B-wf, on the other hand, a steeper spectrum may be attributed to a weaker cascade associated with the lower amplitude (with respect to SIM B-sf) of the injected fluctuations, rather than to dissipation effects associated with ion-heating 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 B-wf, there is no evidence of a clear cascade developing along the direction parallel to B₀ (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.

Fig. 1. For our four simulations, 1D k∥-averaged magnetic energy spectra versus k⊥, normalized so that they overlap at k⊥di ≈ 2.

In Fig. 2 we show for (a) SIM A and (b) SIM B-st, the 3D rendering of the current density magnitude in shaded isocontours; we would like to 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 for the current density J_th, defined in the next section, are shown.

Fig. 2. 3D rendering of regions with J ≥ Jth in blue for (a) SIM A and (b) SIM B-st. 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 current-density 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 $J_{\mathrm{th}}=\sqrt{⟨J^2⟩+3\sigma (J^2)}$, where $\sigma (J^2)=\sqrt{⟨J^4⟩-{(⟨J^2⟩)}^2}$ and ⟨…⟩ denotes the spatial average over the whole simulation box. We want to remark that, by this procedure, any isolated point with an above-threshold 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 grid-point 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:

$$\begin{array}{cc}& \mathrm{elongation}(\mathrm{current}\mathrm{str}.)\mathcal{E}=1-{\ell }_{\mathrm{med}}/{\ell }_{\max }\hfill \end{array}$$(4)

$$\begin{array}{cc}& \mathrm{planarity}(\mathrm{current}\mathrm{str}.)\mathcal{P}=1-{\ell }_{\min }/{\ell }_{\mathrm{med}}.\hfill \end{array}$$(5)

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² (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 $1/\sqrt{{\sigma }_{\min }}$, $1/\sqrt{{\sigma }_{\mathrm{med}}}$, and $1/\sqrt{{\sigma }_{\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, 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 $1/\sqrt{{\sigma }_{\min }}$, $1/\sqrt{{\sigma }_{\mathrm{med}}}$, and $1/\sqrt{{\sigma }_{\max }}$ in the exact same way we did with characteristic lengths and to define two “configuration shape factors” as follows:

$$\begin{array}{cc}& \mathrm{elongation}(\mathrm{magnetic}\mathrm{conf}.)\mathcal{E}=1-\sqrt{{\sigma }_{\min }/{\sigma }_{\mathrm{med}}}\hfill \end{array}$$(6)

$$\begin{array}{cc}& \mathrm{planarity}(\mathrm{magnetic}\mathrm{conf}.)\mathcal{P}=1-\sqrt{{\sigma }_{\mathrm{med}}/{\sigma }_{\max }}.\hfill \end{array}$$(7)

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 point-by-point 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:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{E}}^{\mathrm{HJ}}={[1-\sqrt{h_{\min }/h_{\mathrm{med}}}]}_{\mathrm{cen}}\hfill \\ \hfill & {\mathcal{P}}^{\mathrm{HJ}}={[1-\sqrt{h_{\mathrm{med}}/h_{\max }}]}_{\mathrm{cen}}\hfill \end{array},\end{array}$$(8)

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 “magnetic-configuration shape factors”, ℰ and 𝒫, evaluated through the matrix N at the center of a given current structure:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{E}}^{\mathrm{NB}}={[1-\sqrt{{\sigma }_{\min }/{\sigma }_{\mathrm{med}}}]}_{\mathrm{cen}}={\mathcal{E}}_{\mathrm{cen}}\hfill \\ \hfill & {\mathcal{P}}^{\mathrm{NB}}={[1-\sqrt{{\sigma }_{\mathrm{med}}/{\sigma }_{\max }}]}_{\mathrm{cen}}={\mathcal{P}}_{\mathrm{cen}}\hfill \end{array},\end{array}$$(9)

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:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{E}}^{\mathrm{AV}}={\langle 1-\sqrt{{\sigma }_{\min }/{\sigma }_{\mathrm{med}}}\rangle }_{\mathrm{str}}={\langle \mathcal{E}\rangle }_{\mathrm{str}}\hfill \\ \hfill & {\mathcal{P}}^{\mathrm{AV}}={\langle 1-\sqrt{{\sigma }_{\mathrm{med}}/{\sigma }_{\max }}\rangle }_{\mathrm{str}}={\langle \mathcal{P}\rangle }_{\mathrm{str}}\hfill \end{array},\end{array}$$(10)

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 “knife-blade” 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 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 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.

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 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 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 B-wf, and SIM B-sf (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 B-wf, this discrepancy between the two distributions (viz., generic points versus only points belonging to current over-densities) seems to be particularly emphasized. In the following we give the mean, median, and standard deviation for the corresponding 1D-distributions 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 B-wf, and SIM B-sf 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 high-resolution 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 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 these distributions, obtained from simulations, to the ones obtained using satellite 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 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 “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 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.

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 subset of “generic points” from our simulations. Second, simulation SIM B-sf 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 B-sf 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 B-wf, and SIM B-sf, 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 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 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 B-wf, and SIM B-sf, 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 high-resolution 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 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 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 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 magnetic-field spectrum in Fig. 1). Whether this is actually due to a weak-like versus strong-like 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 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, and standard deviation values in the fifth and eighth rows in Table 3, comparing SIM B-wf and SIM B-sf 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 B-st. As for the previous Sect. 5.1, this very good agreement with SIM B-sf 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 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 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 micro-instabilities) 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 high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019) (see Appendix B for details).

Fig. 9. Normalized occurrence distribution of |b ⋅ emin| (left column) and |j ⋅ emin| (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 tsatA for SIM A, tsatB − wf for SIM B-wf, and tsatB − sf for SIM B-sf. The distribution refers to times at fully developed turbulence, in particular tsatA for SIM A, tsatB − wf for SIM B-wf, and tsatB − sf for SIM B-sf. We superposed these distributions, obtained from simulations, to the ones obtained using satellite 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 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 small-scale 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 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₀ 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₀ 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 sub-group 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⁷. In this case, the subset of “generic” points has been created by considering one point for every 27 (corresponding to a collection of ∼10⁶ 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⁶ (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⁴, 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 tsatA 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 Sorriso-Valvo et al. (2018).

6. Conclusions

In this work we have conducted a wide-spanning analysis of overall nonlocal 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 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 “knife-blade”-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 cigar-like structures (they use the expression “filament-like”) through visual inspection of their 3D Hall-MHD simulations of turbulence. This suggests that ion-kinetic effects (that is, beyond just the Hall term) and/or electron-inertia 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 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 (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 𝒫: “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 𝒫 and elongation ℰ 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), which 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 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 B-wf.

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 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) 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.aida-space.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 Max-Planck/Princeton Center for Plasma Physics (NSF grant PHY-1804048).

References

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 eddy-turnover 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 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, 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.2-0.6, thus they are more filament-like.

Fig. A.1. Scatter plot thickness versus current density for SIM 2D at teddy − 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 ≥ Jth. 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 sub-panel 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 = tsatA 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 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. The choice of these MMS data intervals is motivated by the fact that our simulation setup (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 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 four-spacecraft 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 well-determined 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 inter-spacecraft distance ℓ_SC at each instant and then at ${(\frac{\delta B}{{\ell }_{\mathit{SC}}B})}^2$ setting 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 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 well-determined 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 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 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 inter-spacecraft 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 $max(\frac{3dx}{d_i},\frac{3d\mathit{y}}{d_i},\frac{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 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⊥di ≈ 2.
In the text

	Fig. 2. 3D rendering of regions with J ≥ Jth in blue for (a) SIM A and (b) SIM B-st. 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 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 these distributions, obtained from simulations, to the ones obtained using satellite data (dotted blue line) from the two high-resolution 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 B-wf, and SIM B-sf, 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 high-resolution 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 ⋅ emin| (left column) and |j ⋅ emin| (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 tsatA for SIM A, tsatB − wf for SIM B-wf, and tsatB − sf for SIM B-sf. The distribution refers to times at fully developed turbulence, in particular tsatA for SIM A, tsatB − wf for SIM B-wf, and tsatB − sf for SIM B-sf. We superposed these distributions, obtained from simulations, to the ones obtained using satellite data (dotted blue line) from the two high-resolution magnetosheath intervals analyzed in Stawarz et al. (2019), see Appendix B for details.

In the text

	Fig. 10. For tsatA 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 teddy − 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 ≥ Jth. 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 sub-panel 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 = tsatA for SIM A. Planarity here was computed using the thickness given by the alternative method explained into Appendix A.
In the text