Electron-Only Reconnection in Plasma Turbulence

Hybrid-Vlasov–Maxwell simulations of magnetized plasma turbulence including non-linear electron-inertia effects in a generalized Ohm's law are presented. When fluctuation energy is injected on scales sufficiently close to ion-kinetic scales, the ions efficiently become de-magnetized and electron-scale current sheets largely dominate the distribution of the emerging current structures, in contrast to the usual picture, where a full hierarchy of structure sizes is generally observed. These current sheets are shown to be the sites of electron-only reconnection (e-rec), in which the usual electron exhausts are unaccompanied by ion outflows and which are in qualitative agreement with those recently observed by MMS in the Earth's turbulent magnetosheath, downstream of the bow shock. Some features of the e-rec phenomenology are shown to be consistent with an electron magnetohydrodynamic description. Simulations suggest that this regime of collisionless reconnection may be found in turbulent systems where plasma processes, such as micro-instabilities and/or shocks, overpower the more customary turbulent cascade by directly injecting energy close to the ion-kinetic scales.


INTRODUCTION
Highly accurate in situ satellite measurements of plasma fluctuations and particle distribution functions in the heliosphere are today the primary experimental tool for the study of plasma turbulence (e.g., [1][2][3]) and magnetic reconnection (e.g., [4,5]). In particular, space missions such as Cluster and MMS provide unprecedented opportunities to investigate multi-scale plasma physics, from the ion-kinetic scales down to the electron-kinetic scales, and thereby to constrain theoretical models of kinetic turbulence and reconnection.
Recently, MMS has measured in the Earth's turbulent magnetosheath a series of electronscale reconnection events in which super-Alfvénic electron jets are never accompanied by ion outflows [6]. Those authors dubbed these events "electron-only reconnection" (hereafter, e-rec), to differentiate them from the usual (fast) collisionless reconnection process (e.g., [7][8][9][10]) in which reconnection takes place in an "electron diffusion region" embedded within a larger "ion diffusion region" and in which both species are expelled in a collimated outflow named the "exhaust" (e.g., [11]). While e-rec has been recently shown to occur in kinetic simulations of magnetic reconnection when a single electron-scale CS is ad-hoc initialized [12] or after a quasi-parallel shock [13], it is difficult to explain e-rec in the context of a classical turbulent cascade, in which turbulent energy is transferred conservatively from magneto-fluid scales down to ion-kinetic scales, thereby coupling to the ion dynamics in the usual way [4,[14][15][16][17][18][19][20][21][22][23][24]. This difficulty also holds in scenarios in which the turbulent cascade of magnetic fluctuations is shown to be mediated by the reconnection of ion-scale current sheets (e.g., [20,22,[25][26][27]). Indeed, to the best of our knowledge, simulations in which energy is injected at scales much larger than the ion-kinetic length scales (d i , ρ i being of the same order) show reconnecting CSs with both electron and ion outflows, regardless of the nature of this injection-whether it be continuous forcing or an initial distribution of magnetic and/or velocity fluctuations. This is true also for satellite observations of heliospheric turbulence (now able to resolve the electron scales; e.g., [28]), except for those MMS observations cited above.
In this paper, we use hybrid-Vlasov-Maxwell (HVM) simulations of freely decaying, 2D-3V turbulence to show that e-rec events qualitatively similar to those observed by MMS can develop in an "MHD-scale" turbulent system if magnetic fluctuations are injected sufficiently close to the ion-kinetic scale. These simulations include the Hall and electron-inertia terms in a generalized Ohm's law [29] that captures the decoupling of the magnetic field from the ion dynamics at ion-kinetic scales and allows the complete unfreezing of magnetic flux with respect to the electron fluid motion at the electron-inertial scale. We also demonstrate that simulations of plasmas under the same conditions but with turbulent energy injected farther away from the ion-kinetic scale do not show e-rec, and instead exhibit standard multi-scale reconnection with both ion and electron outflows. The transition from standard reconnection to e-rec is found to occur when the wavenumber range of the injected fluctuations gets close to the ion-kinetic scales, namely when the largest injected wavenumber is varied from (k ⊥ d i ) max = 0.3 to (k ⊥ d i ) max = 0.6. Finally, we show that the physics underlying some features of the e-rec phenomenology can be described by the equations of electron magnetohydrodynamics (EMHD; [30][31][32]).

Basic Equations
We integrate the Vlasov equation for the ion distribution function f i (t, r, v) using an Eulerian approach [33], coupled with Faraday's law of induction for the evolution of the magnetic field B(t, r). In dimensionless units, these are, respectively, Within the hybrid-Vlasov-Maxwell (HVM) framework, these equations are coupled to an Ohm's law for the electric field E(t, r) that includes contributions from the pressure and inertia of an isothermal electron fluid [29]: All equations are normalized with respect to the ion quantities: the mass m i , charge e, inertial length d i , and cyclotron frequency i .
= eB 0 /m i c, where B 0 is the strength of the mean magnetic field. The electron fluid is characterized by a finite electron skin depth d e = √ m e /m i , a flow velocity u e = u i − J/n, and a constant electron temperature T e . The number density n and the ion momentum density nu i are computed as the zeroth and first velocity-space moments of f i , respectively. Quasi-neutrality is assumed, n e ≃ n i = n, and the displacement current is neglected in Ampére's law, so that the current density J = ∇ × B. Electron inertia, see Equation (3), is a key ingredient in this work as it allows for collisionless reconnection to occur even in the framework of a fluid-electron model. At the same time, the Hall term (hidden in the −u e ×B term, equivalent to −u i ×B+J/n×B) is the key actor for speeding up the dynamics around the Xpoint and leading to "fast" magnetic reconnection (see [35] or Appendix B in [34] for further details). We mention that a fully kinetic approach would be more appropriate for the physical description of the diffusive effects since the agyrotropic electron pressure tensor plays a major role (see for instance [36,37]). However, such an approach is out of the scope of this paper and left for future investigation.

Simulation Setup
We consider an initially uniform, Maxwellian, proton-electron plasma embedded in a homogeneous out-of-plane magnetic field, B 0 = B 0êz . The plasma beta parameter for the ion species initially satisfies β i . = 8πnT i0 /B 2 0 = 1, where T i0 is the initial ion temperature. We set T e = T i0 . A reduced mass ratio of m i /m e = 144 is employed, ensuring that d i and d e are well-separated. The HVM Equations (1)-(3) are then solved in a 2D-3V phase space with 1, 024 2 grid points spanning a real-space domain of size L x = L y = 20πd i , corresponding to a uniform spatial resolution ≃0.06d i ≃ 0.7d e . In order to avoid spurious numerical effects at very small scales, we adopt high-order spectral filters [38] on the electromagnetic fields that act only on the high-k part of the spectrum for numerical stability, while reconnection is driven by electron inertia. The velocity space is sampled by 51 3 uniformly distributed grid points spanning [−5, 5]v th,i , where v th,i = √ β i /2 is the initial ion thermal speed, corresponding to a resolution of 0.2v th,i in each velocity direction. We note that velocity fluctuations smaller than the grid resolution, v = 0.2 (in normalized units) in all directions, are well-recovered by the Eulerian approach thanks to the very low level of noise of the algorithm. Indeed, even when the mean flow is small (with respect to v) the asymmetries in the values that the distribution function assumes on the negative and positive part of the vspace grid are well captured by its first-order moment vf (v)d 3 v. We caution that this may not be true only if the width of the distribution is not well resolved, i.e., if v/v th > 1.
We note that velocity fluctuations slower than the grid resolution, dv = 0.2 in all directions, are well recovered by the Eulerian approach where the full sampling of the distribution function allows to capture any mean-flow asymmetry except for the cold plasma limit.
Decaying turbulence is initialized by imposing random, isotropic magnetic-field perturbations, δB = δB ê z + δB ⊥ with ∇ · δB = 0. Such perturbations self-consistently excite "compressive" fluctuations in the velocity and density fields [23]. The choice of injecting compressive magnetic perturbations (rather than, for instance, purely Alfvénic fluctuations having δB = 0) is justified by the fact that inside the bow shock (i.e., in the magnetosheath) the majority of the fluctuations are observed to be compressive (see, e.g., [39,40]). Furthermore, the injection scale in our simulations resides (purposefully) close to d i , at which "Alfvénic" or "magneto-sonic" fluctuations are not clearly distinguishable. Accordingly, we let the (kinetic) response of the plasma to the magnetic-field perturbations to develop velocity and density fluctuations self-consistently. This particular choice has been indeed shown to not affect the development of the turbulent cascade at (and below) the ion scales [23]. However, a study about the precise dependence on a much wider variety of injection properties is beyond the scope of this work-which itself is meant to be a first proof-of-concept-and will be left for future investigation.
Two simulations, labeled sim.A and sim.B for convenience, have been performed. In sim.A, the wavenumbers k of the injected magnetic-field perturbations lie in the range 0.1 ≤ (k ⊥ d i ) inj ≤ 0.6 with a corresponding rms amplitude of the in-plane perturbations |e z × ∇ψ| rms ≃ 0.15. These values are chosen so that moderate-amplitude fluctuations are injected on scales close to the ion-kinetic scales, so as to separate the ion and the electron dynamics from the very beginning of the simulation (see the discussion in section 2.3). In sim.B we excite fewer modes whose wavenumbers lie farther away from the ion-kinetic scales, viz., 0.1 ≤ (k ⊥ d i ) inj ≤ 0.3, corresponding to the regime studied in Franci et al. [20,41], Cerri et al. [23,42], and Cerri and Califano [22]. We adopt the same rms value of the magnetic fluctuations as in sim.A. A third simulation, sim.C, which is equivalent to sim.B in terms of injection scales (viz., 0.1 ≤ (k ⊥ d i ) inj ≤ 0.3) but employs a larger rms value for the fluctuations (viz., |e z × ∇ψ| rms ≃ 0.21), has also been performed. It confirms that the qualitative results from sim.B do not depend strongly on the initial level of the fluctuations, although a more precise statement will require further studies that are left for future investigation. Therefore, since there is no more information in sim.C that is not already present in sim.B, in this paper we only present results from sim.A and sim.B.

The EMHD Limit at Sub-Ion Scales
As turbulent fluctuations cascade to smaller and smaller scales, their characteristic frequencies increase, eventually exceeding those described by the MHD model. At scales below the ionkinetic scales but above the electron skin depth, the ions decouple from the magnetic field, which remains frozen in the electron fluid. In the limit where the spatial derivatives of the number density n and of the ion bulk flow u i can be neglected with respect to the derivatives of the electron fluid velocity u e , the plasma may be described by EMHD. In this case, Equations (2) and (3) reduce to Equation (7) of Bulanov et al. [31], in which the field B ′ .
= B − d 2 e ∇ 2 B (rather than B) is frozen into an incompressible electron flow and the current is carried solely by the electron fluid, so that u e = −J/n. Under these conditions, collisionless EMHD in two dimensions ensures the Lagrangian [31,43], where the flux function ψ is related to the magnetic field via B = e z × ∇ψ + B z . Note that while B ′ and F stay frozen in the electron fluid, B and ψ can reconnect. One advantage of the HVM model adopted here is that it includes EMHD as one of the limits of the generalized Ohm's law (3), while simultaneously accounting also for non-EMHD effects (e.g., compressibility) and the kinetic response of the ions. On this point, we also note that the electron inertia terms in Equation (3), aimed first at allowing reconnection without resistivity, should include a term proportional to d 2 e ∇(∇ ·E). Even if it is negligible in the EMHD limit, this term could influence the kinetic ion response since the electric field directly enters the ions' Vlasov equation. Our set of equations can be seen as the simplest model describing kinetic ions, fluid electrons, ion and electron decoupling at different scales without including dissipation, finally conserving the EMHD flux F in the highfrequency limit. We nevertheless underline that in Equation (3) for the electric field we do not formally impose ∇ · E = 0. As a consequence the electrostatic feedback is anyway included in our hybrid model (see for instance the test on Landau damping on ion acoustic waves in [29]). Numerically, we have two points to solve the electron inertial length d e which is a somewhat borderline to completely separate the EMHD invariant from the flux function. This means that numerical filtering or intrinsic dissipation of the algorithm to advance the Vlasov equation (see [33]), could contribute to the reconnection process and partially dissipate F and so ψ. We note (see section 3.1, in particular Figure 2) that F and ψ are quite well separated during their evolution suggesting that, even at a resolution that must menage at the same time the large-scale evolution and the small-scale reconnection events, electron inertia terms play an important role in allowing for reconnection to occur.

Sim.A: Signatures of e-rec in Kinetic Turbulence
We first describe the results from sim.A. Within roughly 1.5 eddy turnover times of the shortest initial wavelength (∼30 −1 i ), CSs form and reconnection onsets. (Fully developed turbulence is realized only on the eddy turnover time of the outer-scale fluctuations, τ ∼ 100 −1 i .) In Figure 1, top left, we show a fullbox view of the out-of-plane current density J z at t = 42.5 −1 i . For comparison, in the top right frame, we show the same quantity from sim.B (discussed in section 3.2). We observe the formation of several electron-scale CSs, e.g., CS1 (indicated by the rectangle box) at (x, y) = (58, 11), CS2 at (23,11), and CS3 at (21,36). These CSs are characterized by widths of a few d e . More interestingly, their lengths L CS ∼ (2-5)d i are shorter than those typically observed in other turbulence simulations (e.g., [19,20,22,23,44,45]). On these scales, the electron fluid, and thus the magnetic field, is expected to decouple from the ions; it is only deep within these thin CSs that the magnetic field ultimately decouples from the electron fluid. These features are shown in Figure 1, bottom frame, which displays J z (black), the out-of-plane electric field in the frame of the electrons E ′ z . = = E z + (u e × B) z (red line), and E z + (u i × B) z (blue line). (Note that the y-range of the bottom plot for E ′ z has just been shifted with respect to the corresponding plot for sim.A because the current peak in sim.B is negative. However, the length of the y-range shown, |y max − y min |, is the same to facilitate a meaningful comparison). E z + (u e × B) z (red), and the out-of-plane electric field in the frame of the ions E z + (u i × B) z (blue), both vs. y at x ≈ 58d i (corresponding to the vertical dashed line in the accompanying top panel). While the electron fluid is seen to decouple from the magnetic field only deep within CS1 (at y ≈ 11d i , where E ′ z ≈ 0), the ions are approximately decoupled over much of the simulation domain (i.e., E z + (u i × B) z ≈ 0 at all scales). The width of CS1 (either inferred by the condition |J z | > J rms z or by using the more accurate method discussed in section 3.3) is ℓ ≈ 4d e , corresponding to the physical range on which e-rec has been observed to occur in the magnetosheath [6]. Figure 2 presents a zoom-in of CS1 after the onset of e-rec and the formation of an active X-point structure at t = 42.5 −1 i . The top panel reveals a quadrupolar structure of B z near the CS, known to be caused by the Hall term allowing the electrons to locally decouple from the ion neutralizing background and to generate via their in-plane current the quadrupole [11,32,46,47]. Quadrupole signature is today a kind of universal trademark of fast reconnection routinely used not only in simulations but also in space observations to assert the presence of a reconnection event [48,49]. Iso-contours of the flux function ψ (dashed white lines) and of the corresponding EMHD invariant F (black lines) coincide on scales larger than d e , where electron inertia is unimportant. Both quantities are advected by the flow and well conserved everywhere except around the X-point, where ψ locally breaks (and reconnects) on scales ∼d e . This separation of the ψ and F iso-contours is a signature of the EMHD regime (see, e.g., [50]). Indeed, in the bottom panel of Figure 2, the xcomponent of the electron flow u e,x shows the typical electron jet structures coming out from the X-point (highlighted by F, the solid black lines). No evidence of corresponding ion outflows is found in the exhaust around the X-point.
In Figure 3, we show data taken from two virtual spacecraft passing through CS1 along the paths traced by the vertical dashed lines in Figure 2. These trajectories are chosen to be similar to those taken by MMS in Phan et al. [6]. In Figures 3A,F  (given by the condition |J z | > J rms z ). Within these boundaries, there are clear signatures of oppositely oriented super-Alfvénic electron jets, identified as the exhausts (Figures 3B,G), without any noticeable corresponding ion outflows (Figures 3C,H). This feature is in qualitative agreement with MMS measurements [6]. Note that ion jets do not appear even if we move the spacecraft trajectories further away from the X-point location, as would be the case in a "standard" reconnection event. The field-parallel electric field E , E ′ z , and the dissipation rate J · E ′ [51] all depart significantly from zero only across the CS (Figures 3D,E,I,J). Although J · E ′ oscillates, its integral across CS1 is positive, indicating the possibility of a region of net magnetic-to-particle energy conversion (our equations only allow conversion into electron bulk, rather than thermal, energy). All of these features seen in CS1 are observed also in the other CSs.
That the dynamics of these CSs is EMHD-like is supported by Figure 4, which shows spectra of the solenoidal u (sol) and irrotational u (irr) contribution to the in-plane ion and electron velocities and of the in-plane magnetic field during the onset of the first e-rec events (viz., t = 42.5 −1 i ). These spectra demonstrate that the ion flow is almost incompressible in the range d −1 e , while the electron flow remains nearly incompressible across an even larger range. We have also verified that the solenoidal contribution to u e,⊥ around the CSs largely dominates over its irrotational counterpart, and that the "EMHD terms" in our generalized Ohm's law dominate the dynamics (viz., |∇ · (nu i u iz )| ≪ |∇ · (nu e u e,z )| and |nu e,z (∇ · u e )| ≪ |u e · ∇(nu e,z )|; see, e.g., [31]).
Crucially, the properties of the fluctuations shown in Figure 4 at t ∼ 42.5 −1 i do not change character once the fully developed turbulent stage is achieved at t p ∼ 110 −1 i (t p corresponds to the time when the rms out-of-plane current J (rms) z reaches its peak; see, e.g., [44]). After nearly 3 eddy turnover times, corresponding to the fully turbulent regime and by which time the ions could have begun participating in the reconnection dynamics, no sign of ion outflows in sim.A is observed. Figure 5 (left panels) shows a comparison between the out-of-plane and in-plane electron and ion flow velocities at t = 114 −1 i in sim.A; the isocontours of −u e,z ≈ J z /n (top left), u i,z (top right), u e,y (bottom left), and u i,y (bottom right) are shown. All electron quantities exhibit many thin structures at the electron scale, which are well correlated with the CSs traced by J z (not shown here, as it is nearly identical to −u e,z ). On the other hand, all of the ion quantities exhibit much smoother variations and, in particular, are largely uncorrelated with the corresponding electron flows or CSs. As in Figure 2, electron jets are also visible (e.g., in u e,y at (x, y) ≈ (10, 20)d i ; note that here the jets are along y), while no corresponding ion outflows are apparent. Despite these e-rec processes, the fluctuation spectra in the ion-kinetic range of the fully developed turbulence are not significantly affected (when compared to sim.B). The magnetic and electric energy spectra continue to exhibit sub-ion-scale power laws close to −2.8 and −1, respectively, in the range d −1 i k ⊥ d −1 e (not shown here), similar to those found in previous gyro-kinetic, hybrid-kinetic, and fully kinetic simulations (e.g., [19,41,42,[52][53][54][55][56][57][58][59][60][61]) as well as satellite measurements (e.g., [2,[62][63][64]).

Sim.B: Standard Reconnection in Kinetic Turbulence
None of the e-rec characteristics highlighted in section 3.1 are observed in sim.B, in which the magnetic perturbations are injected farther away from d i (at least up to the time when the turbulence is fully developed, t ≈ 200 −1 i ). In particular, as soon as CSs form in sim.B (after 1.5 eddy-turnover times associated with the shortest initial wavelength, ∼60 −1 i here), we observe the development of "standard" reconnection with electron structures embedded in ion macro-layers. This is shown for sim.B in the top-right panel of Figure 1, which provides a fullbox view of the out-of-plane current density J z at t = 95 −1 i . In particular, we see the formation of a CS located at x = 18, y = 50 (highlighted by the black square). As in sim.A, the CS width collapses down to the d e scale at the X-point. However, now its length along the in-plane magnetic field is much larger (more than 10d i ) than those of the CSs observed in sim.A (left panels). As a consequence, at a sufficiently large distance from the X-point, ions couple to the magnetic-field/electron dynamics and ion outflows are generated. This is shown in the bottomright panel (data vs. y taken along the dashed line): by way of comparison with its sim.A counterpart (bottom-left panel), the ions are seen to be better coupled to the magnetic field away from the CS. This facilitates the ions' participation to the reconnection process; indeed, a zoom-in of this CS (see Figure 6) shows the X-point structure with both electron and ion outflows present. Even if smoother, the ion outflows point in opposite directions along the reconnecting in-plane magnetic field and are localized inside the separatrices.
To provide further evidence that the exhausts are composed of two coupled ion-electron jets, in Figure 7 we show six successive plots of the electron and ion outflows measured to the right of the X-point of the CS shown in Figure 6. This sequence shows a central, well-collimated, strong electron jet (continuous lines) that becomes less collimated farther away from the X-point. We also observe a smoother, less intense but clear ion jet superposed on the electron one (dashed lines).
Such an electron-ion structure is never observed in sim.A where the ions decouple from the magnetic-field and the electron dynamics over nearly the entire simulation domain (see Figure 1, bottom left frame, and Figure 5, left frame). This difference is observed not only when the first active CSs form, but also at later times when the turbulence is fully developed. We also note that, although the ion-electron coupling remains valid in the fully developed turbulence stage of sim.B, some electron-scale structures emerge, as shown in Figure 5, right frame, where finer structures develop inside an ion scale reconnecting structure. In summary, in sim.B the different CSs exhibit "standard" or even multi-scale reconnection with both well-developed ion and electron outflows. On the contrary, in sim.A only e-rec develops and "standard" reconnection is never observed, neither at the time when the first CSs activate nor when the turbulence is fully developed.

Statistical Analysis of Current Sheets' Properties
Quantitative differences between sim.A and sim.B can be further assessed by examining the statistical distribution of the CSs' characteristic widths and lengths. For this purpose we define a CS as a region where the magnitude of the current density is bigger than a given threshold J thr . Following Zhdankin et al. [65], we define such a threshold as J thr = < J 2 The widths and lengths of the CSs are then evaluated using two different techniques.
The first technique employs an automated procedure that calculates the width and length of each CS based on the shape of the local |J z |. First, we define the (local) current peak as the maximum value of |J z |. Then, given the eigenvector associated with the largest eigenvalue of the Hessian matrix of |J z | at the peak, we look at the interpolated profile along this direction and define the CS's width as the full width at half maximum of |J z |. Performing a similar procedure to compute the length, namely, interpolating |J z | along the direction of the eigenvector associated with the smallest eigenvalue, could be misleading since the CSs are rarely "straight." For this reason we infer the length as the maximal distance between two points belonging the same CS (i.e., within the thresholded contour), following Zhdankin et al. [65].
We then evaluate the average CS width and length and their standard deviations by taking into account all of the current peaks present in a simulation between the time when the first CSs form (t =  (sim.B). These quantities are shown in Figure 8 for both sim.A (red markers) and sim.B (blue markers): the crosses indicate the average width/length of these runs, while the red (blue) shaded area indicates the parameter region within one standard deviation of these values for sim.A (sim.B). Note that the average width (a few d e ) and its standard deviation are very similar amongst sim.A and sim.B, suggesting that these characteristics are independent of the details of the energy injection. This is not particularly surprising, since the minimum width of a CS width is limited by the decoupling of the magnetic field from the electron dynamics. By contrast, the length distributions are noticeably different in sim.A and sim.B. The average CS length in sim.A is ≃2.6d i and the standard deviation is ≃2.9d i , while in sim.B these values correspond to ≃4.8d i and ≃6.2d i , respectively. Therefore, CSs in sim.B are not only longer on average, but they also show a more pronounced variability in length. The picture emerging from sim.B therefore conforms to the "standard" idea of plasma turbulence developing "a hierarchy of dissipative coherent structures extending down to electron scales" [66], whereas sim.A provides the first numerical evidence of the possibility to "shortcut" this hierarchy and develops only smaller-scale structures where e-rec occurs.
This first analysis procedure has the advantage of being automated, leading to a large statistical ensemble of CSs. That being said, this procedure does not distinguish between reconnecting CSs and non-reconnecting ones. For example, CSs across which the in-plane magnetic field (the only one that can reconnect in 2D) does not change sign but whose current densities are sufficiently large are counted by the automated procedure, despite their inability to reconnect. There are also CSs across which the in-plane magnetic field changes sign, but no clear signs of reconnection (e.g., electron and/or ion outflows emerging from an X point) are present. Our second technique is therefore to compute the widths and lengths of only those CSs that exhibit clear evidence of reconnection. At a fixed time, we visually inspect each CS (identified by our first technique) to assess whether the in-plane magnetic field changes sign and whether electron outflows are present. We then examine those CSs' time evolution to find evidence for the erosion of ψ at the X point, indicating an ongoing topological modification of the magnetic field. The measured widths and lengths of these reconnecting CSs are marked by stars in Figure 8. We have verified, using the proxies J z and E z ′ , that the estimated widths of these CSs represented by a star are correct.
Note that the widths of these reconnection sites are always smaller than the average CSs' width, in both sim.A and sim.B, corroborating the fact that active CSs collapse at the X point. By contrast, the lengths of reconnecting CSs exhibit quite different behavior amongst these two simulations: in sim.A the majority of reconnecting CSs have a length in between the boundaries defined by the standard deviation, with only two cases (out of 8) of reconnecting CSs with a slightly larger length (although smaller than 9d i ). On the other hand, all but one reconnecting CSs in sim.B are out of the boundaries defined by the standard deviation and have a length larger than 10d i . The only one with a smaller length, at ≃2.9d i , is the CS found at t = 203 −1 i where fine electron structures are visible (see Figure 5).
Even if the number of these reconnection sites is statistically smaller than the total number of CSs, we can identify a trend regarding CS width and length within the two simulations. In general, the ensemble of all the CSs in a simulation is not particularly representative of the reconnecting ones. In fact, all the reconnecting CSs have a width that is smaller than the average value of the corresponding distribution of CSs. This width is clustered around ∼ (0.2-0.25)d i ≈ 3d e in both simulations, with a relatively small dispersion around that value (the width being always less than 0.35d i ). More interestingly, the lengths of the reconnecting CSs exhibit a completely different trend amongst the two simulations. In sim.A, the length of these reconnecting structures (which exhibit e-rec) is clustered around ∼ 3.5d i ≈ 42d e , and nearly all of them have a characteristic length which is within one standard deviation from the average value of the distribution of all the CSs in sim.A. In sim.B, on the other hand, nearly all of the reconnecting CSs (which exhibit ion outflows) have a length that is beyond one standard deviation from the average value of the distribution of all the CSs, i.e., larger than 10d i = 120d e .

CONCLUSIONS AND DISCUSSION
Using numerical simulations of kinetic plasma turbulence in a hybrid Vlasov-Maxwell model that includes fully non-linear electron-inertia effects (except for those related to the electron pressure tensor terms), we have demonstrated for the first time that short CSs in which electron-only reconnection takes place can naturally develop within kinetic plasma turbulence at β ∼ 1 when fluctuations are injected on a range of wavenumbers near ion-kinetic scales (in our case, k ⊥ d i = k ⊥ ρ i ≤ 0.6). In this situation, all CSs showing clear evidence of ongoing magnetic reconnection exhibit electron exhausts that are unaccompanied by ion outflows, with the ions being decoupled from the magnetic-field and electron dynamics almost everywhere. The properties of these reconnecting CSs, including their statistical lengths and widths are in qualitative agreement with those "electron-scale" CSs recently measured in the turbulent magnetosheath by MMS [6], in which reconnection-driven electron outflows were unaccompanied by ion outflows. It is worth noticing that in this case where electron reconnection occurs associated with smaller CSs than for standard reconnection (sim.A), we do not observe the classical 2D island coalescence process [67], at least during the full duration of the simulation up to t ≃ 140. In our simplified model in which electrons are modeled as a fluid and decaying turbulence is initialized via compressive magnetic perturbations, this transition from "standard" to "electron-only" dynamics is observed by varying the maximum wave number of the initial fluctuations. In particular, if fluctuations are injected at wavelengths such that k ⊥ d i ≤ 0.3, a hierarchy of dissipative CSs emerge with a variety of lengths, whose evolution follows the "standard" reconnection dynamics in which both ion and electron exhausts are seen.
We have further shown that the wavelength range of energy injection has a measurable impact on the statistical properties of all the CSs: their average length and corresponding standard deviation are indeed different in sim.A and sim.B. In particular, sim.A shows a distribution of CSs' lengths that is centered around ≃2.6d i with a standard deviation of their (asymmetric) distribution of ≃2.9d i , whereas the corresponding distribution in sim.B exhibits a larger average value (≃4.8d i ) and a much wider distribution, with a standard deviation of ≃6.2d i . This difference is even more apparent if we limit the analysis to only those CSs that can be clearly identified as reconnecting ones: in sim.A, these CSs are mostly clustered around a length of ≃3.5d i and are never longer than ≈ 8-9d i . Moreover, their distribution is relatively thin and similar to the one obtained using all the CSs. By contrast, in sim.B we find both very long reconnecting CSs (up to ∼20d i ) and short ones (a few d i ). These statistical properties directly relate to the physics of reconnection that is observed: short CSs in sim.A only produce e-rec, while long and short CSs in sim.B allow for the presence of ion outflows and finer electron structures. Note that, differently to reconnection in a 1D initial Harris sheet [68,69] where the length of the resulting CS is not limited in size and grows in time following the initial configuration, in the turbulent case CSs' characteristic length is limited. Indeed CSs develop in between magnetic flux ropes whose size sets the maximal length CSs can reach. Therefore, the energy injection scale plays an important role by setting the size of flux ropes, and so on the CSs' maximal length and on the possible kind of reconnection developing. Recent fully kinetic simulations of magnetic reconnection developing in an isolated current sheet [12] are in agreement with this picture.
Despite our initial conditions being oversimplified and not directly tied to any specific physical injection mechanism, we are encouraged by the close resemblance of sim.A's results to MMS data. That this resemblance holds not only during the onset of the first reconnection events in sim.B, but also during the fully developed turbulent state, is particularly convincing. We are then led to conjecture that energy injection occurring near ion-kinetic scales, perhaps due to velocity-space instabilities and/or shocks, can qualitatively alter the evolution of CSs and the dynamics of magnetic reconnection in plasma turbulence, potentially explaining the MMS results. In particular, if some energy-injection mechanism with properties similar to those employed in our simulations occurs past the bow shock, our results suggest that e-rec events could be detected relatively close to the bow shock and also further downstream in the magnetosheath. A more recent observational study using an extended MMS data set from the turbulent magnetosheath [70] seems to support this scenario. As the MMS mission goes on, the amount of collected data will enable a more detailed statistical analysis of these e-rec events. This kind of analysis should include also data collected outside of the magnetosheath, in the pristine solar wind, where different statistics for the reconnecting current sheets may be expected. However, due to the different level of fluctuations in the two regions, we stress that a high sensitivity of the spacecraft's instruments may be required in order to detect (the presumably smaller amount of) e-rec events in the pristine solar wind. Next-generation space missions would be most informative.

DATA AVAILABILITY STATEMENT
The simulation dataset (UNIPI e-rec) is available at Cineca on the AIDA-DB. In order to access the meta-information and the link to the raw data, look at the tutorial at http://aida-space.eu/ AIDAdb-iRODS.

AUTHOR CONTRIBUTIONS
All authors have contributed to the data analysis, text editing, and to the discussions presented in the paper.

FUNDING
This project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant