Electron inertia effects in 3D hybrid-kinetic collisionless plasma turbulence

The effects of the electron inertia on the current sheets that are formed out of kinetic turbulence are relevant to understand the importance of coherent structures in turbulence and the nature of turbulence at the dissipation scales. We investigate this problem by carrying out 3D hybrid-kinetic Particle-in-Cell (PIC) simulations of decaying kinetic turbulence with our CHIEF code. The main distinguishing feature of this code is an implementation of the electron inertia without approximations. Our simulation results show that the electron inertia plays an important role in regulating and limiting the largest values of current density in both real and wavenumber Fourier space, in particular near and, unexpectedly, even above electron scales. In addition, the electric field associated to the electron inertia dominates most of the strongest current sheets. The electron inertia is thus important to accurately describe the properties of current sheets formed in turbulence at electron scales.

netic turbulence are relevant to understand the importance of coherent structures in turbulence and the nature of turbulence at the dissipation scales. We investigate this problem by carrying out 3D hybrid-kinetic Particle-in-Cell (PIC) simulations of decaying kinetic turbulence with our CHIEF code. The main distinguishing feature of this code is an implementation of the electron inertia without approximations. Our simulation results show that the electron inertia plays an important role in regulating and limiting the largest values of current density in both real and wavenumber Fourier space, in particular near and, unexpectedly, even above electron scales. In addition, the electric field associated to the electron inertia dominates most of the strongest current sheets. The electron inertia is thus important to accurately describe the properties of current sheets formed in turbulence at electron scales. a) munozp@mps.mpg.de

I. INTRODUCTION
Collisionless turbulence is ubiquitous in space and astrophysical plasmas such as the solar wind, planetary magnetospheres and the interstellar medium. This collisionless (also often called kinetic) turbulence is inherently different in many aspects to the more widely studied and understood fluid or magnetohydrodynamic (MHD) turbulence 1,2 . The differences mainly have to do with the effects associated with the kinetic nature of these plasmas.
Those effects play an important role, in particular, at the dissipation scales of the turbu- as well as high-performance computer simulations, which are capable of resolving a large range of spatial and temporal scales 3,4 .
Turbulence is closely related to magnetic reconnection in different ways 5,6 . Magnetic reconnection is the fundamental process in the Universe by which magnetic energy is quickly transformed into particle energy, often in an explosive way like in solar flares and in the Earth's magnetopause and magnetotail 7 . Magnetic reconnection preferentially occurs in current sheets with a thickness around ion kinetic scales. Since a few decades, the formation of current sheets and magnetic reconnection through them has been observed and analyzed via simulations of turbulence by using a variety of fluid and kinetic plasma models [8][9][10][11][12][13] . This process has also been confirmed by in-situ observations in space plasmas, ranging from the turbulent magnetosheath to the solar wind [14][15][16] .
The very often used inertialess-electron hybrid-kinetic codes, where ions are modeled kinetically while electrons are modeled as a massless fluid, are essential tools to understand plasma turbulence 17 . This is because they can resolve the physics associated to ion-kinetic scales, i.e., the ion skin depth and ion gyroradius in space, and the inverse of the ion cyclotron frequency in time. From the turbulence standpoint, one of the reasons for the importance of those scales is that the nature of the typical collisionless turbulence changes when the energy cascades from scales larger to smaller than the ion-kinetic length scales in space and astrophysical plasmas. This energy cascade is expressed by the spectrum of, e.g., magnetic field fluctuations in the wavenumber (or frequency) space. Above ion-kinetic scales, the turbulence spectrum can be fitted by a power-law with a power-law exponent that can be predicted by MHD theories of turbulence. Below ion-kinetic scales, the powerlaw exponent of the magnetic turbulence spectra suddenly changes, reflecting a change in the turbulent nature of the plasma. This phenomenon is commonly known as "spectral break" in the wavenumber or frequency space. This break has been often observed in-situ in magnetospheric and solar wind plasmas and it often depends on parameters such as the plasma-beta 18,19 . One possible theoretical explanation (among many others) is because at those small scales turbulence is composed by a much larger variety of different plasma waves that depend on ion scales, like kinetic Alfvén waves, in comparison with the pure MHD/Alfvénic-like waves at larger (inertial) scales 20,21 . In addition, hybrid-kinetic models can correctly describe ion-scale current sheets formed in collisionless turbulence, which are typically observed in space plasmas. Because of those facts, among other reasons, there has been a large number of publications about the role of ion-kinetic scales in collisionless turbulence with hybrid-kinetic codes 10,11,13,22 (and references therein).
One important caveat of the usual inertialess-electron hybrid-kinetic approach is that it neglects all electron effects, which are clearly relevant for the correct description of turbulence at electron kinetic scales. There are different kinds of electron effects which can be important for those processes. Some of them, such as the electron Landau damping with its associated dissipation, are purely electron-kinetic effects which can be only captured with a fullykinetic approach for both ion and electrons, or with advanced reduced kinetic models such as gyrokinetics 23 . Other electron effects, potentially relevant for turbulence, are electron temperature anisotropy instabilities, eventually driving turbulence and generating waves 24 .
Temperature anisotropy instabilities result as a consequence of the temperature components evolving differently along different spatial directions, in contrast to a simple but often used in hybrid codes scalar pressure (or temperature). This can be taken into account by means of a hybrid code that implements an evolution equation for the full electron pressure tensor, which in addition also models non-gyrotropic effects (resulting from off-diagonal terms in the pressure tensor) 25 . Effects related to the electron Larmor/gyro-radius are also relevant. They manifest, for example, as finite Larmor radius effects and/or in the form of wave-particle interactions such as electron gyro-resonances which are relevant for maser mechanisms and particle acceleration/diffusion and scattering in turbulence 26 . Finally, the consideration of a finite electron mass (more commonly called electron inertia), is perhaps one of the most important electron effects since it sets a specific length scale, the electron skin depth, which might control another spectral break in collisionless space plasma turbulence at the smallest dissipation scales [27][28][29][30] . As for current sheets and magnetic reconnection through them, electron inertia has three main effects on them: it provides a physical mechanism to break the frozen-in condition in the current sheets (even more efficient under guide-field conditions, see Ricci et al. 31 , Hesse et al. 32 ), it limits the current sheet thinning 33,34 and provides additional pathways to dissipate magnetic energy at the smallest kinetic scales, via, e.g., turbulence generated by electron shear flow instabilities 35,36 . Among those effects, the breaking of the frozen-in condition is perhaps the most important one, since it enables the non-ideal electric field that is necessary for magnetic reconnection and its associated fast magnetic energy release. Current sheets with a thickness around the electron-skin depth have also been recently observed in the turbulent Earth's magnetosheath thanks to the high-resolution measurements of the MMS mission. In those currents sheets magnetic reconnection events without ion-coupling have also been detected 37 .
Although a fully kinetic plasma approach (among other plasma models that include electron kinetics) takes the electron inertia into account, it could be hard to disentangle its effects from the large number of other complex kinetic interactions at electron scales. That is why a hybrid-kinetic plasma approach that takes the electron inertia into account is the ideal tool to investigate those effects separately. Even though some hybrid-kinetic codes have implemented this plasma model, most of them have considered some approximations which might hinder their applicability for scenarios where strong electron-scale gradients are expected, like in turbulence and reconnection.
In the following we briefly review some of those codes which we already discussed in a more extended way in Muñoz et al. 38 and Jain et al. 34 .
The early attempts of hybrid-kinetic codes with electron inertia used the Darwin approximation of the Maxwell equations with the electromagnetic potentials in a Helmholtz decomposition (transverse and longitudinal part) 39,40 . The electron inertia was considered as a term in the transverse (to the magnetic field) electron current density equation and an evolution equation for the electron temperature was utilized. However, they are restricted at least in the sense that the electron inertia is not considered along the parallel (to the mag-netic field) electron equation of motion, and therefore the electron cyclotron frequency is not properly modeled in that direction. Most of the later hybrid codes with electron inertia used a formulation based on electric and magnetic fields. For example, Swift 41  the advantage of a much lower level of numerical noise at the expense of a larger memory requirement in order to store the 6D grid in phase space. The electron inertia was included in a different way, as part of a Helmholtz equation for the generalized electric field that depends on, among other quantities, ion fluids quantities. The latter is a consequence of the use of the CAM (current advance) method used to couple the Maxwell equations with the (ion) Vlasov equation. This involves the evaluation the ion distribution function to get its momenta (ion bulk flow and ion pressure tensor), which are then used to determine the electric field. This code has been mainly applied to solar wind turbulence. Finally, we also mention the code by Amano et al. 44 . The main feature of this hybrid-kinetic code is its capability to model plasma regions with low density and vacuum regions. The electron inertia is included in the generalized Ohm's law for the electric field, but neglecting many of the spatial gradients associated to the currents and densities at ion scales. For a more detailed discussion about the abovementioned codes and others, please see Muñoz et al. 38 and Jain et al. 34 . In particular, a detailed comparison of the different approximations for the electron inertia can be found in Section 2.1 of Muñoz et al. 38 .
We took a different approach to the problem of implementation of electron inertia in a hybrid-kinetic plasma model by developing the hybrid-kinetic PIC (Particle-in-Cell) code CHIEF that implements the electron inertia term of the generalized Ohm's law without approximations 38 . We have efficiently parallelized this code and applied it to the problem of 2D collisionless turbulence, in particular assessing the effects of the electron inertia 34 .
In the present work we extend the work of Jain et al. 34 from 2D to 3D collisionless turbulence. We focus on the effects of the electron inertia on turbulence and the current sheets formed out of it. Many numerical simulations, including our previous study, are performed in 2D configurations in order to save computational resources, even though turbulence is intrinsically 3D. For example, turbulence theories, confirmed by both fluid and kinetic simulations, predict that 2D simulations of this process miss several critical properties compared to their realistic 3D counterparts, such as the dominant waves in the nonlinear energy transfer and modified wave physics that are involved in the anisotropic turbulent energy cascade, resulting in a modified dissipation rate [45][46][47][48] . Current-sheet properties are also expected to be different in 3D compared to their 2D counterparts, as evidenced by the large number of numerical investigations of 3D magnetic reconnection. [49][50][51][52][53] There is thus a clear need for a study of 3D collisionless turbulence with self-generated current sheets and a focus on the electron inertia effects. Our hypothesis is that the turbulence properties should exhibit a change of behavior near the electron skin depth, and that the electron inertia should be relevant for most of the current sheets formed in this turbulence. To the best of our knowledge, this is the first attempt to address this problem with a hybrid-PIC kinetic plasma model that includes electron inertia. Note that there have been previous 3D hybrid-kinetic simulations of turbulence with codes that include electron inertia, but they have not focused on its effects and they were based on a hybrid-Vlasov concluded that the electron Landau damping was responsible for the differences in the spectra between the hybrid and the other two plasma models, and a comparison between the gyrokinetic model with the fully kinetic and hybrid models revealed the importance of fast magnetosonic and Bernstein waves in the turbulent spectra (those two waves modes are ordered out in gyrokinetics). Similarly, in the present article we shall compare two different plasma models in order to reach a conclusion on the contribution of the electron inertia to the current sheets formed out of turbulence.
The remainder of this article is as follows. In Section II A we describe the simulation model. The specific simulation setup is presented in Section II B. Section III shows our results. In particular, in Section III A we show the time evolution of turbulence, useful to identify at which time the analysis are performed. Section III B describes the threedimensional structure of turbulence, so to identify the spatial structure of current sheets.
The core of our new results about electron inertia effects is presented in Section III C.
We discuss our results and state our conclusions in Section IV. In Appendix A we show the parallel scalability of our simulations with the newest version of our code, while in Appendix B we show the variations in the total energy and its components in our simulations.

II. MODEL AND SIMULATION SETUP
A. Hybrid-kinetic plasma model We carry out 3D simulations of decaying turbulence by employing the hybrid-PIC code CHIEF 38 . The hybrid-kinetic model with electron inertia that is implemented in this code has been described in detail in our previous publications 34, 38,60 . Here, we briefly recapitulate the main characteristics of the model: We model ions kinetically via the Particle-in-Cell (PIC) method, as Lagrangian macroparticles whose distribution function f i obeys an ion Vlasov equation. Via an averaging procedure in phase-space we obtain equations of motion for those macroparticles, formally identical to the Newton's equation of motion under the influence of electromagnetic fields, but with the difference that those electromagnetic fields are determined by integrating (interpolating) the full electromagnetic fields by the macroparticles' shape function. The Boris pusher is used to advance those macroparticles in time. The same shape function is used, in turn, to determine the current density (and plasma density) at each grid point from the The electrons are modeled as a fluid employing the quasineutrality approximation n i = n e .
This implies that no continuity equation for the electron density is used, but rather their density is equal to the ion density, which is directly determined from the particle information.
The quasineutrality approximation neglects charge separation and therefore limits our hybrid model to phenomena with a frequency lower than the electron plasma frequency. The electron momentum equation can be rewritten as a generalized Ohm's law for the electric field: Here, the electron inertia effects are contained in the term proportional to m e , which can also be written as − me e due,z dt . ⃗ u e is the electron fluid velocity, so that (⃗ u e × ⃗ B) z is the convective electric field, ⃗ E+(⃗ u e × ⃗ B) the non-ideal electric field. P e,jk is the jk component of the electron pressure tensor. The repeated index k implies a summation, so that ∂P e,jk /∂x k represents the j = x, y, z component. For the sake of simplicity, we used a scalar pressure and a simple isothermal equation of state here so that p e = 1 3 P e,kk = n e k B T e , with k B the Boltzmann constant and T e the electron temperature, constant in time. In Eq. (1), η ⃗ J is the resistive term. We set η, the physical resistivity, equal to zero, but note that there is always some numerical resistivity in the simulations that acts like an additional effective η in the resistive term. This is an effective residual resistivity that arises from the numerical errors of the discretization schemes in the algorithm itself, and is therefore hard to quantify and control.
This way, Eq. (1) shows that the non-ideal electric field in our model can be supported by the resistivity or the electron inertia term. In hybrid-kinetic codes without electron inertia, the non-ideal electric field is only supported by collisional resistivity (either physical or numerical). The mechanisms that support the non-ideal electric field are relevant because they are responsible for breaking the frozen-in condition and allowing magnetic reconnection in current sheets. In collisionless space plasmas the resistivity is not very important, so our model allows for a physical mechanism, the electron inertia, which is appropriate for those conditions and associated to a length scale depending on plasma parameters, the electron inertia.
The electric and magnetic fields obey the Faraday and Ampère's equation without displacement current, which is an appropriate approximation for the considered frequency range of this plasma model: lower than the electron plasma frequency but up to frequencies comparable with the electron cyclotron frequency, thus eliminating light wave propagation. Note that most space plasmas are characterized by an electron plasma frequency which is much larger than the electron cyclotron frequency, making the assumption of both quasineutrality and range of applicability of the hybrid-kinetic plasma model with electron inertia (near and below the electron cyclotron frequency) valid. By combining Eq. (1) with the Faraday's law, we obtain the following evolution equation for the magnetic field in the form of a generalized continuity equation where is the generalized vorticity. Eq. (2) is solved by a parallel flux-corrected transport algorithm.
After the codes solves Eq. (2) for the generalized vorticity − → W , the magnetic field is calculated by combining the definition of − → W in Eq. (3) with Ampère's equation in the following way where ⃗ u i is the ion bulk flow velocity. This is an elliptic equation that is solved by a parallel multigrid method. The electric field is then calculated from the generalized Ohm's law  60 . In particular, the current numerical implementation of this plasma model in our code CHIEF is described in detail in Jain et al. 34 . That publication also describes the parallelization and scaling performance of its (quasi-)2D turbulence simulations. For comparison purposes and to show the scaling behavior of our 3D CHIEF simulations, we discuss those points in the Appendix A.

B. Simulation setup
We use similar initial conditions as those often employed in several previous studies of collisionless plasma turbulence mentioned in Section I. In particular we use a setup similar to the 2D one used by Jain et al. 34 . Note that, following the terminology of This way, the initial magnetic field perturbations, δ ⃗ B ⊥ (with the ⊥ stands for the perpendicular direction to the background magnetic field), are the following randomly-phased, uncorrelated Alfvénic-like fluctuations, where ⃗ k is the wavevector of the fluctuations, ϕ B is a random phase, N kx ,N ky ,N kz are the number of excited modes along each spatial direction, B rms is the root mean square of the magnetic field fluctuation amplitude, chosen as B rms /B 0 = 0.24, which is a value commonly measured in the solar wind 61 . Note that the factor in front of the cosine term in Eq. (5) was chosen in order for each mode to have the same amplitude, i.e, B ⊥ (k x , k y , k z ) = constant.
We excited the longest two wavemodes that can fit in the simulation box and all their is the ion skin depth), and similarly for k y and k z . The ion bulk flow velocities have a similar form as in Eq. (5) but with different random phases ϕ V , so that they are uncorrelated with respect to the magnetic-field perturbations. Their amplitudes are such that the total energy of the ion bulk flow velocity perturbations is the same as the total energy of the magnetic field perturbations, i.e., where V A is the Alfvén speed based on the background magnetic field B 0 . Those initial magnetic field and ion bulk flow velocity perturbations at large (inertial) scales will cascade down to small kinetic scales, eventually forming current sheets.
Other physical parameters are: the background density is constant and equal to n 0 . The electron and ion plasma beta are β e = β i = 2µ 0 n 0 k B T e /B 2 0 = 0.1 with the (scalar) electron temperature T e equal to the ion temperature T i . No initial ion temperature anisotropies were considered. The ion-to-electron mass ratio is m i /m e = 25, which leads to a scale separation between the ion and electron skin depths of d i /d e = 5. The ion gyroradius is The numerical parameters can be summarized as follows. The simulation box spans i . This implies a grid cell size ∆x = ∆y = ∆z = 1.0d e . The timestep is chosen to be ∆t = 0.004 Ω −1 ci . We initially impose a homogeneous distribution of 60 protons per numerical cell. The protons initially have a bulk drift speed given by Eq. (6). We use periodic boundary conditions in all spatial dimensions.
The choice of numerical parameters above (as well as the numerical algorithm itself), in particular the number of particles per cell, can affect the energy conservation of the system, possibly influencing the physical results. One way to improve energy conservation by reducing numerical noise is the use of higher-order shape functions, like the second-order one that was used in our simulations. In order to verify the conservation properties, we show in Appendix B the evolution of the total energy and their components of the simulation Run 1, as well as the calculation method. In general the total energy is conserved relatively well within 1% at the time of the diagnostics. Later there are larger deviations up to 4% towards the end of the simulation, as a result of the dissipation of turbulence.
In addition to the above described main 3D simulation (from now on denoted Run 1 ) we carry out a supplementary simulation in order to analyze the electron inertia. This way our simulations are: • Run 1 : The above-described main 3D hybrid-kinetic simulation with electron inertia (m e ̸ = 0).
• Run 2 : Same as Run 1 but without electron inertia, i.e., the equations that CHIEF solves are the same as those employed by the inertialess-electrons hybrid-kinetic codes (m e = 0).

A. Time evolution
As is customary in collisionless plasma turbulence simulations, a meaningful way to quantify the evolution of the system are the RMS values of the current density. This quantity is associated with the development of current sheets due to the turbulence as it cascades down to ion scales. Fig. 1 shows the relevant component of the current density J rms,∥ ∼ J z (parallel to the background magnetic field) for the standard Run 1. The 3D case without electron inertia (Run 2 ) shows an evolution which is very similar to Run 1 and is thus not displayed here.
J rms,∥ for Run 1 peaks around t ∼ 120Ω −1 ci . At this time the current sheets are quite disrupted, possibly due to reconnection, the formation of secondary structures, and their mutual interactions (see Fig. 2(a3)). A relevant time for analyzing the current sheet structure is therefore before the peak but not too early when current sheets are not fully formed.

C. Electron inertia effects
In order to analyse the effects of the electron inertia term on the current-sheet properties, we first calculate the one-dimensional integrated (reduced) Fourier spectral power of the current density. This is done, in general, by integrating the power of a given scalar field A in the wavenumber Fourier space, A 2 (k x , k y , k z ), on each concentric non-overlapping shell in the perpendicular wavenumber k ⊥ ≈ k xx +k yŷ with a width of ∆k = 2π/L x . This is followed by either an integration along z (parallel direction) in order to get the averaged power as a function of the perpendicular wavenumber, P (k ⊥ ), or an integration along all those shells in k ⊥ -space in order to get the averaged power as a function of the parallel wavenumber,   The results of this calculation for the parallel current density J z (associated with the current sheets) are shown in Figure 4. Although this is not a standard diagnostic for turbulence simulations in comparison with the magnetic field spectra, it helps to assess the scales at which the electron inertia effects influence the current sheets. Note also that the spectral features of the current density J z in Figure 4 do not directly map to those of the magnetic field spectra B ⊥ . So it is not possible to conclude that features (or absence of them) in the spectra J z should also appear in the spectra B ⊥ . For example, spectral power-law indices differ in different regimes, in particular above and below ion-scales (as often investigated in previous works analyzing ion spectral breaks), and the rapid drop at the dissipation scale in the spectra of J z does not appear in the spectra of B ⊥ . Figure 4(a) shows that P (k ⊥ ) is very similar for the two presented cases at ion scales, i.e., for wavenumbers smaller than those corresponding to the electron skin depth kd e = 1 (right gray dashed line). Those spectra are well above the noise level, as shown by the gray curve which represents the initial spectral power of Run 1. Note that the initial bump in the gray curve for low wavenumbers represents the initial Alfvénic wave perturbations. For wavenumbers larger than kd e = 1, P (k ⊥ ) of the current density is larger for the 3D case without electron inertia (Run 2 ) than for the case with electron inertia (Run 1 ) around and below the wavenumber corresponding to the electron skin depth. This implies a flatter spectrum above electron scales for the case without electron inertia. Note that, coincidentally, the differences between Run 1 and Run 2 also appear near the scales associated with the ion gyro-radius, kρ i = 1, since this quantity is close to the electron skin depth for our parameters. However, the differences in the spectra cannot be interpreted as the result of effects related to the ion gyroradius because they only appear when the electron inertia is disabled in Run 2, so they must be related to the latter. Figure 4(b) for P (k ∥ ) shows a similar behavior as for P (k ⊥ ) but with the important difference that the electron inertia effects starts to play a role already at scales much larger than the electron skin depth (approximately at kd i ∼ 2 or kd e ∼ 0.4).
This might imply that the current sheet structure along the parallel direction significantly differs between the cases with and without electron inertia already at scales where electron effects are not normally expected to play a role. Note that this can only be observed in a 3D turbulence geometry; the 2D-limit misses this effect, even though current sheets are also formed.
Note that the power spectrum evolves in time, so it is hard to make precise comparisons for different simulations. Another caveat is that Figure 4 does not allow to pin down whether current sheets are the only responsible structures for the power at electron scales, where differences appear. The following discussion, however, shall provide evidence that this is indeed the case. In order to identify where the main differences caused by the electron inertia appear, we analyze the values of J z in real space. We note that the general excess of spectral power in the parallel current density fluctuations of Run 2 (no inertia) vs Run 1 (inertia), cf. All in all, we conclude that the electron inertia effects tend to limit the maximum values of parallel current density, as well as to increase the volume that is occupied by large values of The previous results indicate that the electron inertia affects the current-density strength, but they do not reveal the causal effects. In order to analyze these, we calculate the contribution of the electron inertia to the electric field in the generalized Ohm's law Eq. (1).
We focus, in particular, on the component of Eq.
(1) that is aligned with the magnetic field, which in our setup is mainly E ∥ ≈ E z :   Fig. 2(a1), we notice that the values of the electron inertia are the largest near some of the strongest current sheets. That statement can be even better justified by the following procedure: First, we identify the strongest current sheets by selecting the grid points with a value of J z > 3J z,RM S , i.e., three standard deviations over the root mean square value. According to Fig. 1, J z,RM S ≈ 0.15J 0 at the time of our diagnostics. The results are shown in Fig. 6(b). Note that this procedure does not allow the identification of all current sheets, since some of them may have current density strengths under the threshold. And this procedure could also pick up values that are not part of current sheets but only strong fluctuations at isolated points due to numerical noise.
Nevertheless, it is a simple and convenient first-order approximation for identifying the strongest current sheets. The second step consists of selecting all the points of the electron inertia term plot in Fig. 6(a) within three standard deviations of the RMS values of the same term (3E RM S,z,inertia ). The results are shown in Fig. 6(c), which highlights the regions where the electron-inertia term is dominant. The third and final step is selecting the values of the electron-inertia term that are located in the regions with strongest current sheets as per the first step, i.e., J z > 3J z,RM S . The results are shown in Fig. 6(d). By comparing with the panel (c) we can thus infer that the electron inertia is dominant in most of the strongest current sheets. But there are also regions with large values of electron inertia that are not located on current sheets, as well as current sheets with values of electron inertia that are not particularly large. Note that this method also has the caveat that it is for a particular 2D plane of the 3D simulation and for a specific time. Analysing data from different 2D planes and different times will lead to current sheets with a different distribution of electron inertia, but the aforementioned trend seems to hold in general. A different way to identify where the electron inertia becomes relevant is illustrated by are not the largest values of current density that are present in the simulation (see Fig. 5).
This value is actually a bit larger than 3J z,RM S . The explanation for this observation relates to the fact that the maximum values of the electron inertia term are located rather at the edges of the current sheets, with corresponding lower values of the current density. This is a known property of the electric field associated with the electric inertia, which has already been found by previous 2D laminar magnetic reconnection simulations 32 . Interestingly, the peak of electron inertia at |J z |/J 0 > 0.5 correlates with the values of current density above which the distribution of currents starts to diverge between Run 1 and Run 2 (3D with and without electron inertia), see Fig. 5. This adds further evidence to our finding that electron inertia plays an important role in sustaining and limiting the growth of most of the strongest current sheets, in particular in a realistic 3D geometry.

IV. DISCUSSION AND CONCLUSIONS
We have carried out 3D collisionless turbulence simulations by using our CHIEF code, which implements a hybrid-kinetic plasma model with a full consideration of the electron inertia term in the generalized Ohm's law. In particular, we assessed the influence of the electron inertia on the current-sheet properties formed out of turbulence. Due to diverse considerations discussed in the introduction, the electron inertia should play an important role on those processes. However, their effects have been so far not investigated in detail in previous numerical investigations of plasma turbulence, even though this phenomenon has been analyzed with plasma models that include those effects, such as fully-kinetic models.
Our previous study has investigated this issue but within a quasi-2D setup 34 . Here we take a step forward by analyzing a more realistic 3D setup. It is known that turbulence is inherently 3D and the 2D limit misses several important aspects, which could potentially affect current sheet properties.
Previous numerical studies of collisionless turbulence have been normally carried out with plasma models that do not include electron inertia, such as MHD or usual hybrid-kinetic models. Those that do include electron inertia, such as the fully kinetic or gyrokinetic plasma models, do not provide an easy way to disentangle the pure electron inertia effects from other electron kinetic effects such as Landau damping. In this sense, the present article consistently compared two different plasma models, hybrid-kinetic with and without electron inertia. This was done within the same numerical code, in order to avoid uncertainties due to the use of different algorithms or numerical schemes.
Our results provide evidence that the effects of electron inertia modify the distribution of parallel current density (associated with current sheets) near electron scales. More specifically, the spectral power of current density is overestimated if electron inertia is not taken into account. This takes place near and below electron scales for the distribution of current density in dependence on the perpendicular wavevector, while at scales larger than electron scales in dependence on the parallel wavevector. This means that the structure of current sheets along the background magnetic field direction is significantly altered if the electron inertia is taken into account. In other words, this modification does not only occurs at electron scales as expected, but also above it, possibly in the range of ion scales. This contradicts our intuitive expectation of electron effects being limited to only electron inertial scales and will be further analyzed elsewhere. Note that this effect can only be observed in a 3D geometry; it is missed in the 2D limit.
Furthermore, we found that the electron inertia also plays an important role in regulating and limiting the largest values of current density. This leads to the formation of current sheets with lower and presumably more physical values compared to plasma models without electron inertia, since a plasma model including electron inertia is more complete than one without it.
Finally, we showed that the electron inertia dominates most of the current sheets. In agreement with previous studies, the maximum values of the electric field associated with the electron inertia are not reached at the peak values of current density, but rather at lower values. Interestingly, those values are the same above which the distribution of currents starts to diverge between the cases with and without electron inertia.
Caveats of our results can be divided in those intrinsic to the model and those related to the parameter choice. We explained both next.
An important caveat intrinsic to the utilized plasma model is the consideration of only the electron inertia as a physical mechanism that can generate a non-ideal electric field and break the frozen-in condition, allowing magnetic reconnection. There are other physical mechanisms that also allow this process; the most important among them is arguably the contribution of the non-gyrotropic electron pressure tensor to the non-ideal electric field. This is neglected in our simulations with the CHIEF code. This mechanism is related to a more accurate description of the electron thermodynamics, which the CHIEF plasma model neglects by assuming a scalar electron temperature governed by an isothermal equation of state. This is a common choice in hybrid-kinetic models and has been shown to reliably reproduce many features of observed space plasmas. Technically, our model could easily adopt, for example, a polytropic equation of state if needed. There are also advanced equations of state used in other hybrid codes aimed at modeling magnetic reconnection, which have proved to reproduce many of the features of this process that can be observed within fully-kinetic plasma models 62 . A more accurate modeling of the electron thermodynamics could use an evolution for the electron pressure (either scalar or tensor), which involves setting additional constants, such as a viscosity 42 . If a tensor pressure is used, its evolution equation also needs to consider additional numerical constants that have to be empirically determined, such as an isotropization term to account for pitch angle scattering due to electron temperature anisotropy instabilities, which are not self-consistently modeled by that evolution equation 63 .
Our choice of parameters regarding the plasma-beta (β e = 0.1) actually also allows to partially justify the usage of the simple scalar equation of state: this low-beta regime, where magnetic pressure dominates over the thermal pressure, implies that the electron gyroradius is smaller than the electron skin depth (ρ e /d e = 0.22). This, in turn, implies that effects related to electron thermodynamics (e.g., equation of state, pressure tensor, associated to the electron gyroradius) should be less important in comparison with the electron inertial effects (which are associated to the electron skin depth). Therefore, a simple scalar equation of state is more appropriate in our case than in higher-beta plasmas.
Regarding other caveats related to numerical parameters, the probably most important one is the choice of a reduced (compared to the real value) proton-to-electron mass ratio.
Our chosen value is m i /m e = 25, which is of course not physically realistic but chosen due to computational resource limitations, as usual in this kind of 3D kinetic plasma simulations.
This could lead to an overestimation of electron inertial effects, however, electron inertia alters the turbulence spectra at electron scales kd e = k 2 ⊥ + k 2 || d e ∼ 1 such that k || << k ⊥ . This is expected to be true for the real mass ratio as well. In any case, our results provide the first evidence for the effects of the electron inertia in this system and represent a first step towards a better understanding of the electron inertia in current sheets formed in kinetic plasma turbulence. Future scaling studies about the dependence of the results on the mass ratio and other parameters have to be carried out in order to establish quantitative predictions which can be directly applicable to space plasmas.
The resolution of our simulations is also not specially high, but sufficient. Indeed, note that the electron gyroradius, as a result of the low plasma beta, is below the grid cell size, ρ e /∆x = 0.22. This quantity, however, does not need to be larger than the grid cell size in our plasma model, which considers electrons as a fluid, and therefore does not include electron gyro-resonances or any related phenomena where the electron gyroradius plays a role. On the other hand, the ion gyroradius needs to be resolved by the grid cell size to properly consider ion kinetic effects such as gyroresonances. In our case, it is resolved by Another different type of concern regarding our results is the level of numerical noise.
This mainly comes from the number of numerical particles per cell, whose value was chosen to be a good compromise between the convergence of values of some physical quantities based on some previous test simulations and the limitations of our computational resources. In order to reduce the numerical noise level, we used a higher-order interpolation scheme/shape function, which smooths out the current density at grid-scales, diminishing the level of numerical noise compared to the more usual linear schemes, and it can contribute to modify the turbulence spectra at wavenumbers corresponding to the grid cell size. In any case, the numerical noise does not affect much our analyzed quantities (in particular current densities) in comparison with other more sensitive parameters to the numerical noise, such as temperatures or electric fields. This is reflected in the low level of numerical noise in comparison with the physical current density spectrum and the good energy conservation properties displayed by our simulations.
In general, we conclude that the electron inertia plays an important tole role to sustain and limit the current density of most of the strongest current sheets formed out of turbulence, in particular around electron scales and more so in a 3D geometry. Unexpectedly, some effects of electron inertia also appear above electron scales, possibly in the range of ion scales. It is thus important to properly take the electron inertia into account in order to better describe the current sheet properties formed out of collisionless turbulence.

ACKNOWLEDGMENTS
The authors acknowledge the support by the German Science Foundation (DFG) projects JA 2680-2-1, MU-4255/1-1, BU 777-16-1 and BU 777-17-1. Computations were performed at the Max Planck Computing and Data Facility (MPCDF). We thank the referees' comments who helped us to improve this paper.

Conflict of Interest
The authors have no conflicts to disclose.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.
In order to assess the parallel scalability of the code, a number of benchmark simulations were performed with our latest version of CHIEF, using the setup of the model described in Section II B. In addition to the setup with N x ×N y ×N z = 256×256×256 grid points (denoted 256 3 for short) used for Run 1 and Run 2, we show, for comparison, also benchmarks for an 8-fold smaller (128 3 ) and an 8-fold larger (512 3 ) number of grid points. In all cases, a number of 60 particles per cell was used, according to Section II B. The performance characteristics of CHIEF for quasi-2D simulations were already documented in Jain et al. 34 .
The computational grid for the 3D simulations and benchmarks is partitioned across MPI processes using a three-dimensional domain decomposition, and assigning an equalas-possible number of grid points and particles to each MPI rank (or processor core). All  For the 256 3 -grid used in this study, the right panel of Fig. 8 shows the strong scaling behaviour of CHIEF with a breakdown into the main algorithmic parts (detailed in Jain et al. 34 , Muñoz et al. 38 ). When running on a single compute node, the overall performance of the code is mainly limited by the available memory bandwidth, which is due to the relatively low algorithmic intensity of 0.77 (a value of at least 20 would be required to achieve peak performance on a single compute node of Raven, according to a basic roofline model) that is achievable by such type of code. When increasing the number of nodes, while the PIC parts of the algorithm scale virtually perfectly in the studied range, the overall strong scalability gets increasingly limited by the less scalable PDE solver in the EMHD part of the code. A detailed analysis has shown that the HYPRE solver is plagued by a significant intrinsic load imbalance beyond a few hundred processor cores and is thus identified as the main scaling bottleneck for the given setup. This is further illustrated by Fig. 9 which shows the run time per time step only for the HYPRE PDE solver as a function of the number of processor cores for different number of grid points, up to a maximum of 640 3 which we could fit into the available memory of the nodes. As expected, the parallel scalability of the HYPRE solver improves (moderately) with increasing number of grid points.
Obviously, shifting the relative computational load towards the well-scaling PIC part by increasing the number of particles per cell will always lead to an improvement of the overall parallel scalability. Therefore, if CHIEF simulations with a number of particles substantially higher than used in this study are desired, the corresponding increase in the (wallclock) run time can be contained by running with a higher parallel efficiency on a larger number (beyond 10, 000) of processor cores. For the simulations of the present manuscript, the total energy is conserved within 1% at the time of diagnostics, while it is within 4% of its initial value at the end of the simulations (See Fig. 10b). The reduction in the total energy is probably due to the dissipation of turbulence. The total energy in Fig. 10 (as well as in our paper with the original description of our code 38 ) is calculated as the sum of the following energy components: where U B is the magnetic field energy, U E is the electric field energy, U i is the total ion kinetic energy, U e is the total electron kinetic energy. Fig.10a) shows how the magnetic energy dominates the other components, obeying the low-beta plasma conditions. The total ion kinetic energy is larger than the electric kinetic energy, while the electric field energy is negligible in comparison with all the other terms. Note in Fig.10b) how the ion kinetic energy is transformed into magnetic field energy and vice versa.
In the following we discuss the calculation of each one of the individual energy compo-nents. The electromagnetic field energies are defined as: where the terms after the first equal sign represent the continuum definition and the terms after the second equal sign represent the discrete approximation used for the calculations.
The indices i, j, k represent the grid points, and N x , N y , N z the number of grid points along each direction in the simulation domain.
In order to determine the kinetic energies U i or U e , we first calculate the ion kinetic energy density per cell u i (the electron energy is obtained by replacing the subscripts i by e) defined as: where p is the particle index running over all the macro-particles in one cell, N p the number of (ion) macro-particles in one cell, and M i is the mass of an ion macro-particle, which is related to the ion mass via M i = (n i ∆x 3 /N p )m i , where m i is the mass of a proton (ion), n i is the physical density, and the factor in front of m i is usually called macro-factor, representing the ratio of physical to numerical particles. The equality from the first to the second row represents a transition from the discrete to the continuum representation. The thermal kinetic energy density and bulk flow kinetic energy densities are then defined from Eq. (B7) respectively as: where the bulk flow energies are defined as and the temperature components as where the index l represents each spatial direction. Note that the terms after the second equal sign in Eq. (B10) and Eq. (B14) represent their discrete version in each cell, but the code calculates those quantities with a weighted average on each cell via a second order shape function.
The total kinetic energy is then the integral of the total kinetic energy density Eq. (B4) over all the simulation domain: n i,j,k k B [T ix,i,j,k + T iy,i,j,k + T iz,i,j,k ] (B14) For our calculations, we used the particle definition Eq. B12 for the total ion kinetic energy shown in Fig. 10. This is approximately equal to the field definition U k,i + U th,i . For