Anisotropic electron heating in turbulence-driven magnetic reconnection in the near-Sun solar wind

We perform a high-resolution two-dimensional fully-kinetic numerical simulation of a turbulent plasma system with observation-driven conditions, in order to investigate the interplay between turbulence, magnetic reconnection, and particle heating from ion to sub-electron scales in the near-Sun solar wind. We find that the power spectra of the turbulent plasma and electromagnetic fluctuations show multiple power-law intervals down to scales smaller than the electron gyroradius. Magnetic reconnection is observed to occur in correspondence of current sheets with a thickness of the order of the electron inertial length, which form and shrink due to interacting ion-scale vortexes. In some cases, both ion and electron outflows are observed (the classic reconnection scenario), while in others -- typically for the shortest current sheets -- only electron jets are presents ("electron-only reconnection"). At the onset of reconnection, the electron temperature starts to increase and a strong parallel temperature anisotropy develops. This suggests that in strong turbulence electron-scale coherent structures may play a significant role for electron heating, as impulsive and localized phenomena such as magnetic reconnection may transfer energy from the electromagnetic fields to particles more efficiently than damping mechanisms related to interactions with wave-like fluctuations.


INTRODUCTION
The solar wind is a highly magnetized and nearly collisionless plasma flow that originates in the Sun's corona (Marsch 2006) and blows in our solar system, filling the the whole heliosphere. While expanding away from the Sun, its temperature gradually decreases with heliocentric distance. Such temperature decrease, however, is smaller than what is expected in the case of a purely adiabatic expansion (e.g., Verscharen et al. 2019, and therefernces therein). This represents the observational evidence that local heating and particle energisation mechanisms are at play in the solar wind (Goldstein et al. 2015). Despite decades-long investigation by different heliospheric missions, understanding the origin of solar wind heating is one of the long-standing open issues in space plasma physics. Among others, turbulent dissipation is the most promising candidate to explain the observed heating (e.g., Kiyani et al. 2015;Parashar et al. 2015).
Turbulence is routinely observed in the solar wind (e.g., Bruno & Carbone 2013). In-situ spacecraft observations return power spectra of magnetic fluctuations exhibiting a power-law behavior over many decades in frequency (e.g., Tu & Marsch 1995;Bruno & Carbone 2013;Chen et al. 2013a;Kiyani et al. 2015). This is the expression of a turbulent energy cascade: the magnetic and kinetic energy stored at the largest scales is transferred (i.e., "cascades") via nonlinear interactions arXiv:2205.08670v1 [physics.space-ph] 18 May 2022 to progressively smaller and smaller scales, down to particles characteristic scales. As the solar wind plasma is collisionless, the energy cascade must terminate via mechanisms that are alternative to particle collisions.
Typically, the power spectrum of the turbulent magnetic fluctuations has a spectral index (i.e., a slope), α, that decreases with increasing frequency. At large scales/low frequencies (the so-called "inertial range"), where the plasma can be described as a magnetized fluid within the framework of magnetohydrodynamics (MHD), both in-situ observations and numerical simulations report α ∈ [−5/3, −3/2], depending on the level of imbalance between outward and inward propagating Alfvén waves and on the heliocentric distance (e.g., Boldyrev et al. 2011;Chen et al. 2013b;Chen 2016;Chen et al. 2020). As we approach particles characteristic scales, where kinetic effects become important, the picture gets increasingly complex. In the range across and just below the ion characteristic scales (known as "spectral break" or "transition region"), both observations and simulations recover α ∈ [−4, −2], depending on plasma parameters such as the ion plasma beta and the level of turbulent fluctuations, α −2.8 being the value more frequently observed (e.g., Bruno et al. 2014;Franci et al. 2016;Verscharen et al. 2019). When the turbulent cascade reaches the electron scales, α further decreases (e.g., Alexandrova et al. 2009;Sahraoui et al. 2013). At present, there is no consensus on the spectral properties of the magnetic fluctuations around and below the electron characteristic scales and even the description of their power spectrum in terms of a spectral index is controversial. In particular, two alternative phenomenological descriptions have been presented, which model the power spectrum of the magnetic field either as a power law with an exponential cut-off (Alexandrova et al. 2009(Alexandrova et al. , 2012Alexandrova et al. 2021), or as two power laws with different slopes at scales larger and smaller than the electron gyroradius (Sahraoui et al. 2013;Huang & Sahraoui 2019). While the spectral index observed in the inertial range can be explained in terms of the Kolmogorov phenomenology (Kolmogorov 1941), eventually extended to include MHD effects (Iroshnikov 1963;Kraichnan 1965;Goldreich & Sridhar 1995), at sub-ion scales we still lack a definite explanation.
In the framework of wave-like turbulence, the energy cascade at kinetic scales is mediated by fluctuations whose properties resemble those of either kinetic Alfvén waves of whistlers waves (Howes et al. 2008;Schekochihin et al. 2009). The mechanisms responsible for dissipating energy and terminating the cascade, thus leading to particle heating, can be either resonant, as for Landau and cyclotron damping (e.g., He et al. 2015;Sulem et al. 2016;Howes et al. 2018;Chen et al. 2019), or non-resonant, as in the case of stochastic heating(e.g., Chandran et al. 2013;Vech et al. 2017;Martinović et al. 2020). Recent numerical studies (Parashar et al. 2010;Papini et al. 2021) have revealed, however, that the turbulent structures at sub-ion scales are characterized by very-low temporal-frequency features, not directly connected to wave activity. This strongly supports the idea that coherent structures such as vortexes, discontinuities, and thin intense current sheets, can play a key role in mediating the turbulent cascade (e.g., Gosling 2007;Greco et al. 2012;Karimabadi et al. 2013;Osman et al. 2014;Franci et al. 2016;Perrone et al. 2016;Wan et al. 2016;Franci et al. 2017;Yang et al. 2017;Camporeale et al. 2018;Papini et al. 2020;Agudelo Rueda et al. 2021;Papini et al. 2021;Smith & Vasquez 2021, and references therein). Particle energization in collisionless plasmas can proceed via the pressure-strain coupling, which starts to play an important role around the ion scales and by which the kinetic energy of a spatially inhomogeneous flow is anisotropically transferred to internal energy (Del Sarto et al. 2016;Yang et al. 2017;Del Sarto & Pegoraro 2017;Matthaeus et al. 2020;Bandyopadhyay et al. 2021;Hellinger et al. 2022).
Another promising mechanism for energy dissipation and particle heating in turbulent plasmas is magnetic reconnection (e.g., Retinò et al. 2007;Gosling 2007;Ergun et al. 2017;Vörös et al. 2017;Eastwood et al. 2018;Stawarz et al. 2019;Zhou et al. 2021). Through reconnection, the energy contained in the magnetic field is released into the plasma, triggered by a topological change in the magnetic field configuration and the subsequent relaxation of the magnetic field lines. Reconnection is usually accompanied by the appearance of ion and electron outflow jets. Recently, however, in-situ observations by Magnetospheric Multiscale (MMS) in the Earth's magnetosheath have revealed some peculiar magnetic reconnection events producing electron jets with no ion counterpart ("electron-only reconnection") ). These were observed in correspondence of thin electron-scale current sheets embedded in strong turbulent fluctuations. Further observational and numerical studies have shown that this is likely to occur when the correlation length of the turbulence is of the order of just a few ion inertial lengths, thus current sheets tend to form and reconnect at scales at which the ions have already decoupled from the magnetic field (Stawarz et al. 2019;Califano et al. 2020;Sharma Pyakurel et al. 2019;Stawarz et al. 2022).
As reconnection converts magnetic energy into particle kinetic and thermal energy, it likely plays a major role in the evolution and dissipation of turbulence (Servidio et al. 2009;Matthaeus & Velli 2011). Indeed, turbulence and reconnection are observed to be intimately linked to each other. Numerical simulations also show that, on the one hand, turbulence develops in reconnection outflows (Karimabadi et al. 2013;Pucci et al. 2017;Lapenta et al. 2020) and the turbulent cascade at sub-ion scales can be directly triggered by reconnection events (Franci et al. 2017). On the other hand, the interaction between turbulent structures spontaneously generates and squeezes thin intense current sheets, until they eventually reconnect (e.g., Matthaeus & Lamkin 1986;Biskamp & Bremer 1994;Wei et al. 2000;Gosling 2007;Servidio et al. 2009Servidio et al. , 2011Franci et al. 2017;Papini et al. 2019;Stawarz et al. 2019;Agudelo Rueda et al. 2021). Magnetic reconnection is also linked to particle temperature anisotropy, as this is expected to be correlated to the hyperbolic flow close to a reconnecting X-point (Cai & Lee 1997;Brackbill 2011;Del Sarto et al. 2016). Anisotropic electron heating and non-gyrotropic velocity distribution functions have indeed been observed both in numerical simulations of single reconnection events (Cai & Lee 1997;Aunai et al. 2013;Pucci et al. 2018a;Sladkov et al. 2021) and in Cluster measurements of a reconnecting thin current sheet in the Earth's magnetotail (Retinò et al. 2008).
Fully-kinetic simulations of collisionless turbulent plasmas, which retain both ion and electron kinetic effects, represent an invaluable tool for investigating the turbulent energy cascade down to electron scales (e.g., Grošelj et al. 2018;González et al. 2019;Roytershteyn et al. 2019;Cerri et al. 2019), its interplay with magnetic reconnection (e.g., Karimabadi et al. 2013;Pucci et al. 2017;Pucci et al. 2018b;Adhikari et al. 2021;Agudelo Rueda et al. 2021), as well as the role of electron-scale coherent structures in dissipating energy and heating particles (e.g., Camporeale & Burgess 2011;Parashar et al. 2015;Yang et al. 2017;Arrò et al. 2021;Bandyopadhyay et al. 2021;Yang et al. 2022). They have also provided numerical evidence for an enhancement of the electron parallel temperature anisotropy in the outflows of strong reconnection events, which occurred spontaneously as the result of the interactions between sub-proton scale turbulent structures (Camporeale & Burgess 2011;Haynes et al. 2014).
In this work, we investigate the development of turbulence, its interplay with magnetic reconnection, and particle heating by means of a high-resolution 2D fullykinetic simulation of strong plasma turbulence in the near-Sun solar wind, carried out using the iPic3D code (Markidis et al. 2010). The initial conditions model the average plasma environment encountered by the Parker Solar Probe (PSP) spacecraft during its first solar encounter at about 36 solar radii from the Sun (Bale et al. 2019). The plasma system size is 32 times the ion inertial length, which allows us to accurately model the turbulent energy cascade from ion scales down to sub-electron scales, providing predictions for the spectral properties of the turbulent fluctuations at scales that cannot be resolved by PSP. To overcome the current unfeasibility of resolving the electron scales while concurrently modelling the large MHD scales, we compare the turbulent spectral properties of the fully-kinetic simulation with those obtained from a hybrid-kinetic simulation (performed with the CAMELIA code Franci et al. 2018), which models a much larger system with a lower resolution, but employing the same physical plasma conditions.
The paper is organized as follows. In Sec.2, we describe the numerical dataset, providing details on the simulation setup and fundamental plasma parameters. In Sec.3, we present our results, in terms of spectral properties of the turbulent energy cascahde, particle heating, and the occurrence of standard and electrononly reconnection events. Finally, we discuss our findings and draw our conclusions in Sec.4.

NUMERICAL DATASET
We performed a 2D particle-in-cell (PIC) simulation of plasma turbulence using the semi-implicit code iPIC3D (Markidis et al. 2010), a fully kinetic code which solves the Vlasov-Maxwell equations for a nonrelativistic plasma of ions and electrons. Ipic3D employs an implicit scheme for the temporal integration of the Vlasov-Maxwell system (Brackbill & Forslund 1982) which removes the numerical stability constraints typical of explicit schemes (Cohen et al. 1989;Hockney & Eastwood 1988), so that it is possible to retain small spatiotemporal scales in an approximate way without the need to resolve the Debye length and to include the speed of light in the Courant condition for the temporal integration. This allows us to retain all kinetic effects, which are vital to describe correctly the overall evolution of the system, while still employing a simulation box whose size is an order of magnitude larger than the ion characteristic scales.
The initial condition consists of a uniform plasma composed of electrons and ions (assumed to be only protons, embedded in and ambient magnetic field B 0 . We use the following normalization units: the magnitude of the ambient field B 0 = |B 0 | for the magnetic fluctuations, the initial plasma density n 0 = n i,0 = n e,0 for the density fluctuations, the (ion) Alfvén speed v A = d i Ω i = B 0 / √ 4πn 0 m i for the velocity fluctuations, the inverse of the proton plasma frequency Ω i = eB 0 /(m i c) for time, the proton inertial length d i = v A /Ω i for lengths. The plasma beta for a given plasma species is β i,e = 8πn i,e K B T i,e /B 2 0 . The electron spatial and temporal characteristic scales are related to the ion ones through the proton-to-electron mass ratio µ = m i /m e as d e = d i / √ µ and Ω e = µΩ i . Quantities and symbols used in these definitions are: the speed of light c, the ion and electron number densities n i,e , the magnitude of the electronic charge e, the proton and electron masses m i,e , the Boltzmann's constant K B , and the proton and electron temperatures T i,e . The 2D computational domain consists of 2048 2 grid points with a spatial resolution ∆x, y = d i /64 and a size L x,y = 32 d i . The time step for the particle advance is ∆t = 0.000625 Ω −1 i = 0.0625 Ω −1 e . We employ a reduced ion-to-proton mass ratio of µ = m i /m e = 100 in order to decrease the computational cost of the simulation and 1024 ions and 8192 electrons per cell to reduce the noise level at small scales. The ambient magnetic field B 0 is along the z-direction. We set its magnitude such that c/v A = ω i /Ω i = 200 and c/v Ae = ω e /Ω e = 20 (we recall here that the electron Alfvén speed is v Ae = B 0 / √ 4πn e m e ). In this work, each field Ψ will be decomposed in its perpendicular (in-plane) component, Ψ ⊥ , and its parallel (out-of-plane, along z) component, Ψ , with respect to the orientation of B 0 . The only exceptions will be the particle temperatures, for which ⊥ and will denote directions with respect to the local magnetic field.
Initially, we assume a uniform number density n = 1, an ion plasma beta β i = 0.2, an electron plasma beta β e = 0.5, and no ion nor electron temperature anisotropy, i.e. A i,e = T (i,e)⊥ /T (i,e) = 1. We add an initial spectrum of in-plane Alfvénic-like magnetic and ion bulk velocity fluctuations, composed of modes with wavenumbers in the range −0.8 k x,y d i 0.8 and random phases. Their initial global amplitude can be estimated in terms of their root-mean-square value as δB rms 0.39 B 0 and δu rms i 0.34 v A , so that there is a small initial residual energy (excess of magnetic over ion kinetic energy). Initial electron bulk velocity fluctuations with both perpendicular and parallel components are also imposed, such that n 0 (δu i, − δu e, ) = −n 0 δu e, = J = ∇ × δB ⊥ and n 0 (δu i,⊥ − δu e,⊥ ) = J ⊥ = ∇ ⊥ × δB = 0. As a consequence, at t = 0 we have δu e,⊥ = δu i,⊥ , so that δu rms e,⊥ 0.34 v A .

Fully developed turbulence at sub-ion scales
We let the plasma system evolve from its initial condition until a turbulent energy cascade has fully devel- oped and starts to slowly decay. This occurs when the rms value of the current density, J , reaches a maximum/plateau (e.g. Franci et al. 2015aFranci et al. , 2017. Figure 1 shows the time evolution of some fundamental spaceaveraged quantities that allow us to follow the development of the turbulent cascade, identifying different phases. Figure 1(a) shows the rms of the components of J both parallel and perpendicular to the ambient magnetic field B 0 . The former provides the dominant contribution to the total current and reaches a maximum at t peak = 11.4 Ω −1 i , marked by a green dashed vertical line. The following analysis will therefore be performed at this time. The other component, J rms ⊥ , is initially almost zero and quickly grows, yet remain-ing quite smaller than its parallel counterpart at all times. This is related to the rapid development of parallel magnetic fluctuations, as shown below. The eddy turnover time corresponding to the initial injection scale, k inj d i (k 2 x + k 2 y ) d i 1.1 is the time associated to nonlinear energy transfers at t = 0 and can be estimated 3 Ω −1 i . A turbulent energy cascade develops fully from the injection scale down to the electron scales in a few nonlinear times, as t peak ∼ 5 t NL . This is similar to what previously observed in hybrid simulations with similar initial conditions (e.g. Franci et al. 2015b,a), where about 10 nonlinear times were required. Figure1(b) shows the time evolution of the maximum of the magnitude of the current density computed over the whole simulation box. This starts increasing after about one t NL and then reaches a first maximum at t rec 4.9 Ω −1 i , which is after about 2 t NL . In Franci et al. (2017), we have provided numerical evidence of the link between the time of the first maximum of |J | and the onset of magnetic reconnection events. Here we observe the same, as we see clear signs of reconnection in the magnetic and current structures starting at t t rec (not shown). Figure 1(c) shows the rms of the in-plane and out-of-plane components of the magnetic field (red solid and dashed, respectively) and of the ion bulk velocity (light blue). The behaviour of B rms ⊥ and u rms i,⊥ is qualitatively the same as previously observed in hybrid simulations (cf. Fig. 1 of Franci et al. (2015a)): after a quick re-adjustment of the initial conditions, the former slightly increases before slowly decreasing, while the latter decreases during the whole evolution. As a result, the initial small excess of magnetic energy further increases reaching a maximum at about half the simulation and then remains almost constant. The parallel components δB rms and u rms i, , which are initially zero, quickly start to increase until they reach an almost constant and comparable value, which is a few times smaller than their perpendicular counterparts. This indicates that the levels of compressibility and magnetic compressibility that spontaneously form are relatively small but not negligible. It is interesting to note that δB rms ⊥ reaches a maximum at t = 5 Ω −1 i t rec and starts decreasing when reconnection events start occurring. Figure 2 shows the isocontours of the energy in the magnetic fluctuations, |δB| 2 = |B−B 0 | 2 = δB 2 x +δB 2 y + δB 2 z (panels (a) and (b)) and in the current density, |J | 2 , both in the whole 32 d i × 32 d i simulation domain (panels (c) and (d)) and in a 8 d i × 8 d i sub-domain (bottom panels). Coherent magnetic structures, i.e., vortexes or islands, are embedded in a more chaotic environment where stretched and twisted shapes emerge. The vor-texes have radii from a few times d i down to a few times d e . Gradients of the magnetic field occur at smaller scales, as the strongest current structures have a width of the order of d e . For an easy comparison by eye, we have also drawn circles with radius r = d i and r = d e in fig. 2(a)-(c) and lines with length l = d i and l = d e in fig. 2

Spectral properties of the turbulent fluctuations
Figure 2 has shown that the turbulent structures are randomly oriented, as expected given that the ambient magnetic field is orthogonal to the simulation plane so there is no privileged direction in the plane. We can then reasonably assume the two-dimensional spectra of the turbulent fluctuations to be statistically isotropic and analyse the spectral properties of each field by computing its omnidirectional power spectrum P. The top panel of Fig. 3(a) shows the power spectra of the fluctuations of the magnetic field B, of the electric field E, and of the ion and electron bulk velocities u i and u e at t = t peak . Black vertical dashed lines mark the wavenumbers corresponding to the ion and electron inertial length, d i,e , and gyroradius, ρ i,e . For each field, a dashed line with the same color indicate the noise level, estimated by computing the power spectrum at t = 0 when only the large-scale initial fluctuations should be present. The bottom panel of Fig. 3(a) shows the signalto-noise ratio (SNR) for each field. In the following, we will consider the power spectra reliable (meaning that their shape is not affected by the noise) at those scales where SNR > 3. This holds for k ⊥ d i 20 (or, equivalently, k ⊥ d e 2) for P E and P ui and for k ⊥ d i 45 (k ⊥ d e 4.5) for for P B and P ue . This ensures that our modelling of the magnetic field power spectrum is meaningful around the electron scales and for a further half a decade in wavenumber above k ⊥ d e = 1. The top panel of Figure 3(a) suggests that both P B and P ue may exhibit a power-law behaviour with different spectral indices in different ranges above and below the electron scales. On the contrary, P ui starts dropping at k ⊥ d i 1 and its behavior is not so clear. There is a hint that for 4 k ⊥ d i 10 it may exhibit a similar slope to the magnetic field (although such interval is too limited to make any claim) and then it starts flattening just before the noise level is reached, so it is not possible to infer whether this is a numerical or physical effect. Finally, P E at sub-ion scales is much less steep than P B , so that the level of the electric field fluctuations is much larger than the magnetic field's for k ⊥ d i 5. In order to provide a more quantitative characterization of P B and P ue , Fig. 3 shows that the magnetic field power spectrum can be modelled by P B ∝ k α B ⊥ with α B −11/3 for almost a decade at k ⊥ d e < 1 and −5 for a shorter interval at k ⊥ d e > 1. Correspondingly, we observe P ue ∝ k αe ⊥ with α e −5/3 and −3 in the two intervals of scales, respectively. The relation α e = α B + 2 can be easily explained by considering the definition of the current density, J = n i u i − n e u e , and Ampère's law which links it to magnetic field, J = ∇×B (where the displacement current has been neglected). Due to charge neutrality, n i = n e = n and, assuming that n = n 0 + δn with δn n 0 , we get u e ∝ ∇×B. In Fourier space this reads u e,k ∝ k × B k , from which we obtain P(u e ) ∝ k 2 ⊥ P(B) (since in 2D k = k ⊥ ). At smaller scales, k ⊥ d e 2.5, there is a hint of another possible power-law range in P(u e ), withα e −5. This interval is very narrow, just about a factor of 2 in wavenumber, as the power spectrum quickly reaches the noise level. Although we cannot consider this power-law behavior and its slope reliable enough, the further steepening of the power spectrum is evident from the top panel of Fig. 3(a). In the same interval, the magnetic field power spectrum also steepens accordingly, although α e = α B + 2 does not seem to hold perfectly anymore, possibly due to the effect of the large noise in the density due to which the assumption δn n 0 might not hold anymore and/or to the fact that at large wavenumbers (and frequencies) we cannot neglect the displacement current anymore.
The fact that at the smallest scales P(u e ) seems to exhibit still a power-law behavior P(u e ) while P(B) seems to follow approximately, might represent a hint that the turbulent cascade is still proceeding at k ⊥ ρ e 1, but the dynamics might be dominated by the electron processes rather than by the magnetic fluctuations. Investigating this will require new simulations with an even higher resolution and larger number of electrons per cell, in order to further improve the accuracy of the power spectra below the electron scales. As mentioned in Sec. 1, there is no consensus in the literature on the shape of the magnetic field power spectrum across the electron characteristic scales and, more precisely, on whether it exhibits a double power-law behavior or rather an exponential decay. So far, we have discussed how two power laws with a spectral index of −11/3 and −5 seems to provide a good modelling for P B over a little more than a full decade across the electron scales, more specifically for 0.25 k ⊥ d e 4. We also compensated the magnetic field power spectrum by the inverse of an exponential cut-off similar to that suggested by Alexandrova et al. (2021), which reads P B ∝ k −8/3 ⊥ exp (−c k ⊥ ρ e ) with c = 1.8. The result is shown by the green line in the top panel of Fig. 3(b). This also seems to provide a good approximation for P B as it is almost a horizontal line, although with significant oscillations, for a little more than a decade in wavenum-ber. We need to note, however, that we have set c = 1 instead of c = 1.8 as in Alexandrova et al. (2021), since the latter does not work well here. The difference in this constant might be due either to the unrealistic mass ratio employed in the simulation or to the fact that the exponential cut-off model fails in our case, possibly due to the specific plasma conditions. Providing strong numerical evidence in favour or against this goes beyond the scope of this paper, as we only intended to show that the simulation is able to model the turbulent cascade accurately up to scales smaller than the electron gyroradius, regardless of what the exact shape of the power spectra is.
A limitation of this fully kinetic simulation, as for most simulations employing the same model, is the fact that the MHD scales are not retained. Usually, especially when explicit numerical codes are used, this is due to the fact that all characteristic spatial scales need to be resolved, including the Debye length. This is not a requirement in our case, as the iPic3D code employs the implicit moment method. Here, however, the box size is still limited as we decided to prioritize employing a very large number of ppc to reduce the numerical noise. While this choice assures that the power spectra are more accurate at large wavenumbers, missing completely the MHD range makes the turbulent cascade start developing at the ion scales. Some important properties (e.g., the plasma and magnetic compressibility, the residual energy) at ion scales are set by our initial condi- . Direct comparison between the power spectra of the kinetic simulation, P kin , and those of the hybrid simulation P hyb for the magnetic field B (a), the electric field E (b), the ion bulk velocity ui (c), the electron bulk velocity ue (d), the ion density ni(e), and the parallel magnetic field B (f). For each field, the ratio P hyb /P kin is also shown in the bottom panels.
tions rather than self-developing from the larger scales. One could then wonder whether this sub-ion scale turbulence is representative of the plasma dynamics in a larger plasma system or if it is instead sensibly constrained or affected by the box size. In order to investigate this and validate our results, we have run a hybrid simulation with the same physical conditions but a much larger box, able to model the turbulent cascade over two decades in wavenumber, one above and one below the ion scales. This simulation has the same parameters and initial conditions as the one described in (Franci et al. 2020a), so that further information can be found in Appendix A therein. The only (minor) differences are that the grid points are 4000 × 4000 instead of 4096 × 4096 (and correspondingly the box size is 250 d i instead of 256) and 2048 ppc instead of 1024.
In Fig. 4 we compare the power spectra of different fields for the fully kinetic simulation, P kin , and for the hybrid simulation, P hyb . For the former, we only draw the portion of the power spectra where the corresponding SNR is above 3. For the latter, since we have only the initial power spectra for u i and n i , we apply the SNR criterium by using the noise threshold SNR(u i ) < 3 for P ui and SNR(n i ) < 3 for all other power spectra. Fig. 4(a) compares the spectral behavior of the magnetic fluctuations. P hyb B is more extended at large scales, since the hybrid simulation box is larger. This allows for the development of a Kolmogorov-like power law with a spectral index of −5/3 for k ⊥ d i 3.5, which then steepens to a power law with a slope compatible with −11/3, as already observed in Franci et al. (2020a). P kin B , instead, completely lacks the −5/3 range and is higher at the largest scales. This is likely due to two main reasons. First, the injection scale is just in the middle of that range, and probably the fluctuations at the largest scales did not have time to fully partake in the cascade. Secondly, in the fully kinetic simulation the mechanisms leading to a fully developed inertial range are missing and therefore those underlying the sub-ion scale turbulence may start determining the dynamics at slightly larger scales. Between the ion-scale break of the hybrid simulation (the end of MHD inertial range) and the electron scales, i.e., in the range 3.5 k ⊥ d i 10, P hyb B and P kin B overlap. Finally, at k ⊥ d i 10, i.e., k ⊥ d e ∼ k ⊥ ρ e ∼ 1, when P kin B steepens as the electron physics kicks in, P hyb B maintains the same slope until the noise level is reached. In Fig. 4(b), P hyb E and P kin E are compared. These are quite overlapped in the range 0.6 k ⊥ d i 6 while, at smaller scales, the former flattens while the latter keeps going down with a more or less constant slope. This difference in behavior may be partially due to the different level of noise in the two sim-ulations, as ∆x hyb = 4 ∆x kin and ppc hyb = 1/4 ppc kin . The latter condition is related to the fact that in the hybrid case the electrons are treated as an isothermal fluid, so that P e ∝ β e /2 n i and the number of ions per cell (2048) determines the noise in the electric field, while in fully kinetic case P e ∝ n e T e , so the number of electrons per cell (8192) is the determining one. As a consequence, a sub-domain of a given size contains 16 times more particles in the fully kinetic simulation with respect to the hybrid one. It is, however, reasonable to ascribe the flattening at scales much larger than the spatial resolution to different physics, given that in the hybrid model we approximate ∇ · P e ∼ β e /2∇n i , so any possible effect due to electron temperature anisotropy and nongyrotropy, which can be expected to become important when the electron scales are approached, is not retained in the hybrid model. In Fig. 4(c) we compare the power spectra of the ion bulk velocity. These are very close to each other at all the scales where the noise is negligible, with just a minor difference when reaching the electron scales. From this comparison, we can conclude that the spectral behavior of the ion velocity seems to be the same in the two simulations, thus it is not significantly affected by the presence of kinetic electrons and related processes. The power spectra of the electron bulk velocity, shown in Fig. 4(d), start overlapping at k ⊥ d i 0.5 and they proceed together for about a decade, up to k ⊥ d e 0.5, where P hyb ue further steepens, accordingly with P hyb B as explained above. Figure 4(e) compares the power spectra of the ion density fluctuations. We note here that for the fully kinetic simulation we have verified that the spectrum of the ion density, P ni , and that of the electron density, P ne , are exactly identical up to k ⊥ d i 40, where they start being affected by the two different levels of noise (due to the different number of pcc for the two species). We observe that P hyb ni and P kin ni almost perfectly overlap in the whole range k ⊥ d i 1 and k ⊥ ρ e 1, with P hyb ni getting larger just below the electron scales. Finally, Fig. 4(f) compares the power spectra of the parallel magnetic fluctuations. Interestingly, P hyb Bz is slightly larger than P kin Bz at all scales, especially below the ion-scale break where they exhibit the same spectral slope therefore keeping a constant differing factor. This could be somehow related to the different initial amplitude of the turbulent fluctuations with respect to the ambient field, given that the ratio P hyb Bz /P hyb Bz is compatible with the ratio B hyb rms /B kin rms 1.2. Another possible reason for the difference in the magnetic compressibility could be related to the plasma betas: although we start the two simulations with the same values of β i and β e , these evolve differently and reach different values when the turbulence is fully developed (in particular, β e does not change in the hybrid case by definition). This difference, however, is very small and as such cannot be consider as an indication of a different regime or different underlying physics. Fig. 4 globally returns a picture where the MHD-scale fluctuations and their properties and evolution do not seem to have significant effects on those below the ion scales. This is supported in particular mode by the spectral behaviour of the density fluctuations, since these are zero at the beginning (apart from the small-scale noise due to the finite number of ppc) and only develop self-consistently through plasma processes. This further validates the findings by Franci et al. (2020b) on the kinetic turbulent cascade: the plasma dynamics in the kinetic range is mainly controlled by a few fundamental physical parameters, so that the properties of the fluctuations at sub-ion scales are independent of the presence and extension of an inertial range.
Since in the fully kinetic simulation β e = 0.5, then ρ e = √ β e d i 0.7 d i so the two electron characteristic scales are very close to each other. It is not possible to infer whether one of the two, or a combination of them, determines the steepening of the power spectra around the electron scales. Preliminary analysis on another fully kinetic simulation (not shown here) with β e = 0.04, where ρ e = 0.2d e , seems to suggest that the electron-scale break is related to k ⊥ ρ e rather than to k ⊥ d e . At this level, however, we cannot exclude that the break might occur at either of the two scales depending on the value of β e , just like it is observed to happen for the role of β i for the ion-scale break (Chen et al. 2014;Franci et al. 2016).

Anisotropic ion and electron heating
We now focus on the electron and proton heating generated by turbulence. In Fig. 5, we report the overall evolution of both ion and electron temperatures and their spatial distribution once turbulence is fully developed. Fig. 5(a) shows the time evolution of the average value (over the whole simulation box) of the perpendicular, T i,⊥ , parallel, T i, , and total temperature, T i = (2 T i,⊥ +T i, )/3, and of the temperature anisotropy, A i = T i,⊥ /T i, , all normalized to their initial values. We recall here that the temperature components are defined with respect to the local mean magnetic field. At the time of maximum turbulent activity, T i has increased of almost 8% with respect to its initial value T in i = 0.1. Such increase does not stop at t peak , being linear in time afterwards. This hints at the possibility that, in the presence of an energy injection mechanism that keeps feeding the turbulent cascade and maintains it in a quasi-steady state over a longer time, the ion temperature could increase significantly more. T i,⊥ starts increasing as soon as the simulation begins and keeps increasing throughout the whole evolution, with some oscillations until t 5 and then a linear increase in time. T i, , instead, remains constant until t 2 (comparable to the eddy turnover time at the injection scale), then decreases reaching a minimum just before t rec and it increases monotonically afterwards. Correspondingly, A i increases until t rec and then it decreases very slowly, reaching an almost constant value of 1.06 at t peak . The corresponding quantities for the electrons are shown in Fig. 5(b). T e remains almost constant until about t rec , then increases very slowly, reaching a value of 1.01 T in e at t peak , with T in e = 0.25 being its initial value. Despite the quite small increase in the total electron temperature, its perpendicular and parallel components exhibit very large variations. T e,⊥ shows a small initial decrease, then stays almost constant until just before t rec and then starts decreasing more rapidly, linearly in time. On the contrary, T e, increases initially more or less linearly in time, then at around t rec it starts increasing much faster. Correspondingly, A e decreases monotonically throughout the simulation (faster after t rec ), reaching a value of about 0.9 at t peak , which keeps decreasing afterwards.
Summarizing, Fig. 5(a)-(b) show that at t t rec = 4.9 Ω −1 i (i.e., when the first maximum of |J | marks the onset of magnetic reconnection), there is a clear and significant change of behavior in the particle heating, characterized by the fact that: (i) the rate of variation of the total ion temperature increases and becomes almost constant; (ii) the total electron temperature starts increasing linearly in time; (iii) the ion temperature anisotropy reaches a plateau and starts decreasing very slowly; (iv) the rate of variation of the electron temperature anisotropy increases. It is then reasonable to assume that the change in particle heating is related to the onset of magnetic reconnection events.
In order to test this assumption, Fig. 5(c)-(h) show the spatial distribution of the above quantities in the 2D simulation domain. In the top row (panels c-e), we show the ion perpendicular and parallel temperatures and their ratio. T i,⊥ and T i, exhibit a very patchy behaviour, with alternating regions where they are larger or smaller of their initial value. Typically, in regions where one component has increased the other one has decreased. As a result, despite their average values are very similar and their maximum values are also comparable, A i is larger than 2 in many large areas, reaching a maximum as large as 3.7. Regions with values well above or below 1 are alternating and somewhere we observe a quadrupolar configuration, both in cor- respondence of X-points where magnetic lines reconnect (e.g., at the point [8,17]) and inside big vortexes (e.g., at [19,9]). The spatial distribution of the electron temperature looks quite different. T e,⊥ is close to its initial value in most of the simulation box, exhibiting just a small decrease in correspondence of what look like reconnection outflows around X-points and a small increase at the center of vortexes. On the contrary, T e, exhibits much larger variations, slightly decreasing in regions where the magnetic field lines get compressed and becoming about twice larger in the reconnection outflows. As a conse-quence, in these latter areas A e gets very small reaching values as small as 0.36. The spatial distribution of both ion and electron temperatures is consistent with what we have inferred from the time evolution of their average values and with the observed change of behavior at the time when reconnection starts occurring: reconnection is likely to be at least partially responsible for particle heating, especially for the electrons. This sounds reasonable, as magnetic reconnection is a mechanism that transfers energy from the magnetic field to the particles.

Standard and electron-only reconnection
In order to confirm the role of reconnection in heating the particles and in generating a strong temperature anisotropy, especially for the electrons, we now provide a more quantitative characterization of the regions where reconnection events occur. Figure 2 clearly showed qualitative evidence for the presence of reconnection events in the simulation, e.g., X-points and small magnetic islands emerging in proximity of thin strong current sheets. Magnetic reconnection is known to occur in turbulent plasmas as the the result of the interaction of turbulent eddies which leads to the formation and stretching of intense current sheets. Indeed, it has been observed to occur spontaneously in numerical simulations of plasma turbulence performed with different methods (e.g. Matthaeus & Lamkin 1986;Servidio et al. 2009;Franci et al. 2015b;Wan et al. 2015;Haggerty et al. 2017) and it has been shown to provide a significant contribution to the further development of the turbulent cascade at sub-ion scales (e.g. Franci et al. 2017;Mallet et al. 2017;Papini et al. 2019). In Fig. 6(a), we show a pseudocolor plot of the out-ofplane component of the vector potential, A z , with its isocontours superimposed as black lines. Green stars mark saddle points, i.e., X-points which represent potential reconnection sites. Fig. 6(b) shows arrows marking the direction of the local magnetic field in the simulation box, on top of the contour plot of the |δB| 2 . This clearly shows that the local magnetic field has opposite sign at the two sides of most of the X-points, as required for magnetic reconnection to occur. Recently, Agudelo Rueda et al. (2021) have investigated the spontaneous occurrence of reconnection in a 3D fully kinetic simulation of plasma turbulence by means of a set of indicators based on thresholds. Here we chose to apply the same criteria to detect and characterize reconnection events. Although our 2D setting is definitely less realistic than a 3D one, here we have the advantage that it is straightforward to identify X-points in 2D. In this sense, our analysis of reconnection events can be considered as complementary to the work by Agudelo Rueda et al. (2021), and in some sense a further validation. For a given (positive) scalar quantity, Ψ, we define a threshold Ψ th = Ψ + N Ψ rms and we look for regions where Ψ/Ψ th − 1 > 0. Here, we chose to set N = 1, as we find that this is enough for detecting the strongest reconnection events in our simulation. The criteria we have evaluated are: C1 : |J | > |J | th ≡ |J | + |J | rms (Fig. 6c) C2 : |u i,e | > |u i,e | th ≡ |u i,e | + |u i,e | rms (Fig. 6d, e) C3 : T i,e > T th i,e ≡ T i,e + T rms i,e (Fig. 6g, h) C4 : J · E| + > J · E| th + ≡ J · E| + + J · E| rms + (Fig. 6f) where for the last criterion, we consider only positive values of J · E, which denote energy transfer from the fields to the particles. Figure 6(c) shows the result of indicator C1 for |J |, revealing the presence of intense current structures. We can clearly see that all these have a width of the order of d e and lengths of the order of a few times d i . We count about a dozen of them, taking into account that some of them are actually portions of the same current sheet, as the simulation box is periodic. In Fig. 6(d)-(e), we estimate the indicator C2 for |u i,e |, related to the presence of fast ions and electrons that provide evidence for the presence of a reconnection outflow (or jet). The ion bulk velocity exceeds the threshold in some regions with a width of the order of 1 − 2 d i , some of which seem to be directly related to X-points and represent ion jets in the reconnection outflow, while others seem to be within large-scale vortexes. On the contrary, the electron bulk velocity is above its threshold in all the regions directly connected to the strong current sheets observed in panel (c). These fast electron jets have a very small width, of the order of the d e , as they are confined in small regions between larger-scale structures. It is interesting to note that many of the observed reconnection events seem to be highly asymmetric, as jets are observed only on one side of the X-point and not on the other. Figure 6(f), shows the contour plot of indicator C4. J · E| + is above its threshold in a few thin regions corresponding to the most intense current sheets, providing further evidence that in the location of reconnection events magnetic energy is converted into particle energy. In Fig. 6(g)-(h), we estimate indicator C3 on T i,e which is related to the presence of heated ions and electrons. The regions where the particle temperatures exceed their respective threshold look similar for ions and electrons, meaning that most reconnection events are eventually leading to both ion and electron heating. A major difference, however, can be appreciated by comparing with Fig.5(c)-(h): reconnection heats the electrons preferentially in the parallel direction with respect to the local magnetic field and the ions in the perpendicular direction. The latter result appears to be an agreement with the fact that a sheared in-plane ion flow tends to transfer energy to the in-plane ion pressure tensor components via the action of the symmetric part of the strain tensor (Del Sarto et al. 2016;Yang et al. 2017;Del Sarto & Pegoraro 2017;Matthaeus et al. 2020;Bandyopadhyay et al. 2021;Hellinger et al. 2022). Again, we observe a significant lack of symmetry, as particles are heated only on one side of the X-point, typically the same for both species. In Fig.6(i), we compare the indicators C2 for ions and electrons, by showing the contour plots of their perpendicular components, |u i,⊥ | 2 and |u e,⊥ | 2 on top of each other (in shades of blue and red, respectively). This allows us to appreciate if a reconnection event produces both in-plane ion and electron jets, as in "standard" reconnection, or only an electron jet with no ion counterpart, as in "electron-only" reconnection. For some of the X-points, it is not straightforward to apply this classification, mainly because no strong ion nor electron jets are observed or because the geometry is complex. We still observe about a dozen of events which can be labelled as electron-only (yellow circles) and about half a dozen that are standard (light blue circles). Summarizing Fig. 6, we can conclude that as a result of the interaction between turbulent magnetic structures, thin intense current sheets form and shrink until they undergo reconnection. Such reconnection events transfer magnetic energy to the particles, heating the electrons in the direction parallel to the mean field and the ions mainly in the orthogonal direction. Thin electrons jets are also observed in the reconnection outflow, with ion counterparts only in about one third of events. The fact that standard reconnection is observed only in a minority of events might be related the simulation box size, which is only a few tens of d i , and to the energy injection, which occurs at the ion scales (Califano et al. 2020). These conditions strongly limit the magnetic correlation length, which characterizes the size of the interacting turbulent magnetic structures and, as consequence, also set the length and thickness of the current sheets forming in-between (e.g. Stawarz et al. 2019). In other words, there are only a few vortexes with radius of the order of few times d i , so there are only a few chances to develop strong current sheets with such a length. On the contrary, most of them are formed by the interaction of magnetic structures with smaller size and will therefore be shorter. Indeed, the ion jets observed in Fig. 6(i) seem to be concentrated just next to some of the largest vortexes. Quantitatively characterizing the properties of reconnection and its interplay with electron-scale turbulence goes beyond the scope of this work and will be the subject of future investigation. Here, we intended to show evidence for the coexistence of both types of reconnection events spontaneously driven by turbulence at sub-ion scales and at electron scales.

DISCUSSION AND CONCLUSIONS
In this work, we have presented the results of a twodimensional simulation of plasma turbulence at sub-ion and electrons scales performed with the fully kinetic code iPIC3D, showing signature of turbulence-driven magnetic reconnection and related anisotropic particle heating. We have modelled plasma conditions similar to those measured by PSP during its first solar encounter. This allowed us to extend towards larger wavenumbers the analysis of the ion-scale turbulent cascade performed with hybrid simulations in Franci et al. (2020a). Our simulation employs a large box in terms of the electron characteristic scales, i.e., 320 d e × 320 d e (d i = 10 d e , since the proton-to-electron mass ratio has been set to 100), a spatial resolution ∆x = d i /64 0.16 d e , and implements a very large number of particles (1024 ppc ions and 8192 ppc electrons, for a total of more than 40 billion particles in the whole simulation domain). This setting allowed us to accurately model the development of the turbulent cascade and to determine the spectral properties of the magnetic and electron bulk velocity fluctuations for a decade and a half in wavenumber in the range 0.1 k ⊥ d e 4, fully capturing the electronscale transition.
The turbulent dynamics is characterized by the concurrent presence of coherent magnetic field structures in the form of vortexes, with radii between a few times d e and a few times d i , and very elongated thick filaments, with a width of the order of d i and a length up to about about 20 d i . In between these, thin intense current sheets form, with width of the order of d e .
Our results show that, at the maximum of turbulent activity, the power spectrum of the magnetic fluctuations exhibits is very well modelled by a double power law, with a spectral index α B compatible with −11/3 and −5 respectively above and below the electron characteristic scales. Since the electron beta is close to 1, the electron inertial length and gyroradius are very close to each other, thus making it impossible to infer whether the electron-scale transition is associated to either or both of them. Complementary numerical simulations (not shown here) performed by varying both spatial resolution and number of ppc, have showed that the location of such transition is of physical origin, provided that a sufficient number of grid points allows to cover approximately a full decade across the electron scales. The power spectrum of the electron bulk velocity behaves just as one would expect considering that at subion scales the current density is almost exclusively due to the electron bulk motion: it also behaves like a double power law, with a spectral index α e = α B + 2. We observe a hint of a further steepening in the power spec-trum of the electron bulk velocity fluctuations (and correspondingly, of the magnetic field) at k ⊥ d e 2. However, such steepening occurs very close to the scale at which the spectra reach the noise level, therefore it is not possible to draw any conclusion. Simulations with a better spatial resolution will be needed to further investigate whether such steepening is present or not. The spectrum of the electric fluctuations behaves differently, as it does not exhibit any major change of slope when the electron scales are reached. Such behavior seems to be consistent with Cluster observations in the Earth's magnetosheath (Matteini et al. 2017), where the steepening of the spectrum of the magnetic field at electron scales does not have a counterpart in the electric field. A thorough investigation of the nature and of the spectral properties of the electric field, including an analysis of the different contributions to the generalized Ohm's law, is currently ongoing and will be the subject of a follow-up work.
The limited size of the simulation box in terms of d i , together with the fact that the initial energy injection in the simulation occurs at the ion scales, makes it impossible to model the turbulent cascade in the inertial MHD range. We have compared the power spectra of all the fields with those obtained from a hybrid kinetic simulation (performed with the CAMELIA code) which implemented the same plasma conditions (same ion beta, same eletron beta, very similar initial amplitude of the turbulent fluctuations) but used a much larger box. Indeed, all power spectra from the fully-kinetic simulation are well in agreement with their hybrid counterparts, down to scales comparable or slightly larger than d e . Below this scale, almost all the spectra from the iPic3D simulation further steepen, while the ones from CAMELIA exhibit no change of behavior. This is not surprising, as the two models differ significantly at the electron scales: CAMELIA treats the electrons as an isothermal massless fluid and therefore electron inertia effects and electrons kinetic effects are not retained. What is not necessarily obvious is the fact that the spectral behavior below d i is almost completely unaffected by the existence/absence of a turbulence cascade in the inertial range. This hints at a certain degree of universality of the sub-ion scale turbulent cascade: its properties seem to depend on the plasma conditions alone, regardless of what is happening at larger scales. This is in agreement with what has been recently observed in Franci et al. (2020b). There, the spectral properties of turbulence driven by a Kelvin-Helmholtz instability observed by MMS in the Earth's magnetosheath have been correctly modelled by a hybrid simulation of Alfvénic turbulence. This employed the observed values for the plasma beta, the electron beta, and the level of turbulent fluctuations with respect to the ambient magnetic field.
The development of turbulence is observed to have an effect on particle heating. The ions quickly develop a certain temperature anisotropy, whose average value reaches a maximum of 1.08 at t t rec and then decreases very slightly. Their average parallel temperature with respect to the local magnetic field starts increasing at t t rec and reaches the average perpendicular component at t t peak . This behavior is quite similar to what observed in previous hybrid simulations, where the perpendicular ion temperature quickly starts to increase, leading to the formation of temperature anisotropy, with a comparable increase of the parallel temperature starting later (more or less at time of the the onset of magnetic reconnection), therefore freezing the temperature anisotropy (Franci et al. 2015b(Franci et al. , 2017. The electrons develop a temperature anisotropy in the opposite sense, i.e., the parallel component is much larger than the perpendicular and whose average value keeps decreasing throughout the whole evolution, even after the turbulent cascade has fully developed. Both the ions and the electrons are heated during the simulation, as their total temperature increases of 7% and 1%, respectively. In the time evolution of the temperature of both species we observe a link with the onset of magnetic reconnection: at t t rec , the rate of increase of the ion temperature sharply increases, while the electron temperature starts increasing after remaining at its initial value in the first part of the simulation. The spatial distribution of the electron temperature and its anisotropy confirms a major role of magnetic reconnection in heating the electrons: the regions where the parallel electron temperature increases (and therefore where the temperature anisotropy mostly differs from 1) are localized in the outflows of reconnection events at the sides of X-points. This has been confirmed by detecting reconnection events using the criteria defined and applied in Agudelo Rueda et al. (2021) (there for a 3D fully kinetic simulation), which are based on intensity threshold on the current density, the ion and electron bulk velocity and temperature, and the energy transfer from the electromagnetic fields to the particles mediated by J · E. The results of such analysis have revealed the coexistence of "standard" and "electron-only" reconnection events, exhibiting both ion and electron jets or electron jets alone, respectively. There is a hint that the different nature of reconnection events is related to the length of the reconnecting current sheet, which in turns depends on the size of the interacting magnetic vortexes (Stawarz et al. 2019). The latter is also directly link to the correlation length of the turbulence, which has been observed to have an impact on the nature of the reconnection (Califano et al. 2020;Sharma Pyakurel et al. 2019;Stawarz et al. 2022). The areas where we observe the largest increases in electron temperature are all located in reconnection outflows. For the ions, some areas of increased temperature are in correspondence of reconnection events, while others appear to be inside magnetic field vortexes.
In conclusion, we have modelled the development of plasma turbulence accurately from ion down to sub-electron scales, observing the spontaneous occurrence of magnetic reconnection events (both standard and electron-only), which are associated with strongly anisotropic electron heating.
In order to make the simulation computationally feasible, in this work we employ smaller values of some plasma parameters with respect to their real value or to their typical value in the solar wind at 1 AU: m i /m e = 100 instead of 1836 (as a consequence, d i = 10 d e instead of 43 d e ), c/v Ai = ω i /Ω i = 200 instead of ∼ 10 3 and c/v Ae = ω e /Ω e = 20 instead of ∼ 200. This is consistent with the fact that here we are modelling the near-Sun solar wind, as encountered by PSP during its first perihelion. In this sense, our parameters are just intermediate between the typical solar-wind values mentioned above and those in the upper solar corona, i.e., c/v Ai ∼ 10 2 and c/v Ae ∼ 3. Verscharen et al. (2020) demonstrated that plasma models employing m i /m e = 100 and c/v Ai 10 can successfully cover physics on scales 0.2 d i for β i ∼ β e ∼ 1, which is the regime we have explored here. Since we employ that exact value of m i /m e , an order of magnitude larger value of c/v Ai , and we also have c/v Ae 10, we are confident that our fully kinetic simulation provide quite a correct modelling of the plasma dynamics down to electron scales.
That said, a main limitation of our study is the 2D geometry, which strongly constrains our simulation results, since both turbulence and magnetic reconnection are inherently 3D processes. Modelling the spectral behavior reliably from ion to sub-electron scales and estimating particle heating quantitatively, however, require a high accuracy (in terms of both grid size and number of ppc) which are difficult to reach in 3D at the moment. Therefore, high-resolution 2D fully kinetic simulations represent an optimal starting point for this kind of analysis. Still, all the results obtained in the present study will need to be validated by future 3D simulations with similar plasma conditions. Indeed, this will be the subject of future work.