Issue 
A&A
Volume 653, September 2021



Article Number  A156  
Number of page(s)  16  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202140279  
Published online  28 September 2021 
Bridging hybrid and fullkinetic models with Landaufluid electrons
I. 2D magnetic reconnection
^{1}
Physics Department “E. Fermi”, University of Pisa, Largo Bruno Pontecorvo 3, 56127 Pisa, Italy
email: francesco.finelli@phd.unipi.it
^{2}
Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA
^{3}
Centre for Mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium
^{4}
Istituto per la Scienza e Teconologia dei Plasmi, Consiglio Nazionale delle Ricerche (ISTPCNR), Via Amendola 122/D, 70126 Bari, Italy
^{5}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire J.L. Lagrange, Bd. de l’Observatoire, C.S., 34229, 06304 Nice, France
Received:
2
January
2021
Accepted:
5
July
2021
Context. Magnetic reconnection plays a fundamental role in plasma dynamics under many different conditions, from space and astrophysical environments to laboratory devices. Highresolution in situ measurements from space missions allow naturally occurring reconnection processes to be studied in great detail. Alongside direct measurements, numerical simulations play a key role in the investigation of the fundamental physics underlying magnetic reconnection, also providing a testing ground for current models and theory. The choice of an adequate plasma model to be employed in numerical simulations, while also compromising with computational cost, is crucial for efficiently addressing the problem under study.
Aims. We consider a new plasma model that includes a refined electron response within the “hybridkinetic framework” (fully kinetic protons and fluid electrons). The extent to which this new model can reproduce a fullkinetic description of 2D reconnection, with particular focus on its robustness during the nonlinear stage, is evaluated.
Methods. We perform 2D simulations of magnetic reconnection with moderate guide field by means of three different plasma models: (i) a hybridVlasovMaxwell model with isotropic, isothermal electrons, (ii) a hybridVlasovLandaufluid (HVLF) model where an anisotropic electron fluid is equipped with a Landaufluid closure, and (iii) a fullkinetic model.
Results. When compared to the fullkinetic case, the HVLF model effectively reproduces the main features of magnetic reconnection, as well as several aspects of the associated electron microphysics and its feedback onto proton dynamics. This includes the global evolution of magnetic reconnection and the local physics occurring within the socalled electrondiffusion region, as well as the evolution of species’ pressure anisotropy. In particular, anisotropydriven instabilities (such as firehose, mirror, and cyclotron instabilities) play a relevant role in regulating electrons’ anisotropy during the nonlinear stage of magnetic reconnection. As expected, the HVLF model captures all these features, except for the electroncyclotron instability.
Key words: plasmas / magnetic reconnection / instabilities / methods: numerical
© ESO 2021
1. Introduction
Weakly collisional, magnetized plasma dynamics are characterized by a wide range of scales, both in space and time. Moreover, the nonlinearities unavoidably at play in collisionless plasma dynamics can simultaneously couple multiple scales within the system. On top of that, plasma processes taking place at kinetic scales are able to provide important feedback on the global (i.e., fluidscale) evolution of a system, for example in the context of magnetic reconnection (MR). Reconnection is probably the best example of a nonideal, microscale process able to change the global magneticfield topology of a system by locally breaking the magnetic connections ideally preserved at magnetohydrodynamic (MHD) scales (Biskamp 1994). The change in magneticfield connections occurring in relatively short times (to the extent that it is sometimes defined as “explosive”) in turn triggers the release of a huge amount of energy stored in the previous magnetic configuration (e.g., Treumann & Baumjohann 2013, and references therein). As for MR, several other processes exist that can lead to the coupling between disparate spatiotemporal scales, from fluid to kinetic. Amongst others, one of the most natural examples is plasma turbulence, which could be defined as a (fascinating) “multiscale disorder” (Schekochihin et al. 2008). These two processes, turbulence and MR, are indeed intimately related to each other, one process possibly feeding on the other and vice versa (e.g., Servidio et al. 2011; Zhdankin et al. 2013; Cerri & Califano 2017; Franci et al. 2017; Pucci et al. 2017; Comisso & Sironi 2018; Muñoz & Büchner 2018, and references therein). In general, the intrinsically nonlinear, multiscale nature of collisionless space plasma dynamics requires a numerical approach to the solution of these problems. This highlights the need for a compromise between the complexity of an appropriate theoretical framework and the available computational resources.
Numerical studies are based first on the choice of the plasma model used to describe, at best, the plasma dynamics of interest while keeping the computational cost of a simulation under control. The fullkinetic approach, based on the solution of the Vlasov equation for all plasma species (viz., ions, and electrons), represents the best mathematical model for the study of the dynamics of a collisionless, magnetized plasma at kinetic scale. However, fullkinetic simulations usually come at a very large computational cost, too large to afford realistic plasma parameters and/or to include a wide enough range of spatiotemporal scales. On the opposite limit in terms of completeness of the plasma description and, as a consequence, of the computational cost of their numerical solution, one has “fluid” models. In this context, the idealMHD theory represents the simplest case of an infinitely conducting, singlefluid plasma where all kinetic effects at any characteristic scale are ignored. In between these two extremes, with respect to the full Vlasov approach, there is a long list of hybrid models (Tronci & Camporeale 2015). The first hybrid models were proposed a very long time ago to address quasineutral plasma dynamics, and they rely on a kinetic description of the ions and a fluid description of the electrons (e.g., Byers et al. 1978; Hewett & Nielson 1978; Winske 1985). The hybrid approach has garnered a great deal of attention over the past two decades because it represents a very good compromise between the physics to be described and the computational cost of the numerical simulations. The capability of hybrid models to catch the correct physics across the ion kinetic scales and below, approaching the electron scales, depends on the model adopted for the electron response. The simplest case of massless isothermal electrons may be too simplified when applied to cases in which the electron response may have important feedback on the ion dynamics. In this respect, within the hybrid approach, it is possible to improve the electron response by first including finite electron inertia effects in Ohm’s law (e.g., Valentini et al. 2007; Muñoz et al. 2018). A step in this direction was taken in Le et al. (2016), where a hybrid model with a static and gyrotropic equation of state for a massless electron fluid was considered. However, that equation of state was derived considering only the adiabatic trapping of electrons as the primary source of anisotropy. A further step can be taken by introducing a more general fluid model that includes dynamic electronpressure anisotropy (Chew et al. 1956; Hunana et al. 2019b), electron finiteLarmorradius (FLR) effects, and/or electron Landau damping (e.g., Hammett & Perkins 1990; Sulem & Passot 2015; Hunana et al. 2019a).
In general, hybridkinetic models have proved to be capable of satisfactorily catching the main kinetic physics at play for a large number of problems, ranging from fluid and kinetic instabilities (Hellinger & Matsumoto 2000; Matteini et al. 2006; Califano et al. 2008; Henri et al. 2013; Kunz et al. 2014) to collisionless shocks (Lembège et al. 2009; Caprioli & Spitkovsky 2013; Weidl et al. 2016), dynamo effects (Rincon et al. 2016; StOnge & Kunz 2018), MR (Le et al. 2016; Palmroth et al. 2017; Cerri & Califano 2017; Franci et al. 2017; Wang et al. 2019; Califano et al. 2020), and kineticscale turbulence (Servidio et al. 2015; Grošelj et al. 2017; Cerri et al. 2018, 2019; Hellinger et al. 2019; Wang et al. 2019). The main goal of our project is to investigate the possibility of improving the electron description in hybridkinetic models. The starting point is the “hybrid Vlasov–Maxwell” (HVM) code (Mangeney et al. 2002; Valentini et al. 2007), which is already equipped with a fluid model in which electrons are described as an isotropic, isothermal fluid with finite inertia. The HVM code has recently been upgraded in order to implement a more sophisticated model for the electron fluid. This includes evolution equations for the anisotropic (gyrotropic) electron pressures, p_{∥, e} and p_{⊥, e} (where ∥ and ⊥ refer to the local magneticfield direction, b = B/B), and a Landaufluid (LF) closure for the parallel transport of the gyrotropic electron thermal energy along field lines (i.e., parallel heat fluxes, q_{∥, e} and q_{⊥, e}). Hereafter, we refer to this new model as “hybrid Vlasov–Landaufluid” (HVLF). The idea is to include within a hybrid description the relevant electron pressureanisotropy effects and a fluid model for the electronkinetic response that still holds in a nonlinear regime. Therefore, the LF model implemented in the HVLF code goes beyond the early attempts to include these effects in simplified settings (e.g., Hammett & Perkins 1990; Snyder et al. 1997; Passot & Sulem 2007) and is based on the approach presented by Sulem & Passot (2015).
In this work we consider the problem of MR, focusing on its nonlinear evolution. Since MR is a typical process in which both ion and electronkinetic effects at subion scales drive an energy conversion that can feedback into the “large” (fluid) scale of the system, it represents a suitable problem in which to employ the HVLF model. In fact, as highlighted by both simulations and observations, the heating associated with this process is typically anisotropic, exhibiting an enhancement of the parallel heating near the region where electrons are demagnetized, the socalled electron diffusion region (EDR; see, e.g., Chen et al. 2008; Daughton & Karimabadi 2007). Magnetic reconnection naturally drives pressure anisotropies, which in turn constitute a freeenergy source for secondary instabilities, most notably firehose instabilities (FHIs) and mirror instabilities (MIs). Their development bounds the plasma distribution in parameter space during the nonlinear stage of MR, as also seems to be the case in the turbulent solar wind (SW; e.g., Hellinger et al. 2006; Matteini et al. 2007; Bale et al. 2009). Here an MR simulation performed with the HVLF model is compared with equivalent simulations (i.e., identical initial setups) performed with (i) the standard HVM code with isothermal electrons and finite electroninertia effects and (ii) a fully kinetic model employing the semiimplicit particleincell “iPic3D” code (Markidis & Lapenta 2010; Lapenta 2012).
In Sect. 2 we introduce the three different models that have been used in this work and the initial simulation setup. Numerical results are then detailed in Sect. 3, starting with an analysis of the reconnection dynamics in Sect. 3.1 (e.g., reconnected flux and reconnection rate, as well as the properties of the relevant fields near the EDR). In Sect. 3.2 we focus on the properties (and possible role) of the pressure anisotropy that is developed by protons and electrons during the nonlinear stage of MR. In Sect. 4 we summarize the results and outline our conclusions.
2. Kinetic plasma models and MR setup
In this section, we present the three different kinetic models that have been used in this work. These are (i) a hybridkinetic model, where the electrons are approximated either (a) as isotropic and isothermal fluid or (b) as anisotropic fluid with an LF closure, and (ii) a fullkinetic model. It should be noted that the hybridkinetic approximation also keeps finiteelectronmass effects in both cases.
It is important to point out that, while classic hybridkinetic models exclude electron kinetic effects, such as electron Landau damping (see, e.g., Told et al. 2016; Camporeale & Burgess 2017), this downside is partially cured when an LF type of closure is adopted for the electron fluid (e.g., Hammett & Perkins 1990; Sulem & Passot 2015).
2.1. Hybridkinetic models
The hybridkinetic approximations consists of a quasineutral plasma, n_{p} = n_{e} ≡ n, where fully kinetic ions (protons in this case) are coupled to a fluidelectron model (Winske 1985). The evolution of the proton distribution function, f_{p}(x, v, t), follows the Vlasov equation,
while the magnetic field, B, evolves following Faraday’s law of induction,
where e and m_{p} are the proton electric charge and mass, respectively, while c is the speed of light. The proton dynamics is coupled to the electron fluid through a generalized Ohm’s law that provides the electric field, E,
where d_{e} is the electron inertial length, the corresponding terms in (3) taking into account finiteelectronmass effects, while u_{p} and u_{e} ≡ u_{p} − J/en are the proton and the electron fluid velocities, respectively. The current density J is computed from the magnetic field via the Ampère’s law, 4πJ/c = ∇ × B, where the displacement current has been neglected in the hybrid approximation, while n and u_{p} are obtained as velocity moments of the proton distribution function.
The set of Eqs. (1)–(3) needs to be closed by one (or more) equation(s) for the pressure tensor of the electrons, Π_{e}. Two different closures for the electron fluid are considered here, isothermal or LF electrons, and the corresponding equations are detailed in the following subsections.
For the numerical solution of (1)–(3), we adopt an upgraded version of the Eulerian HVM code (Valentini et al. 2007), which now allows different electron closures. When a standard isothermal electron fluid is considered, we refer to it simply as HVM (Sect. 2.1.1), while when the new LF electron model is adopted, we refer to it as HVLF (Sect. 2.1.2). In the HVM code, an Eulerian approach is employed to numerically solve the Vlasov equation (1) for the proton distribution function (i.e., f_{p}(x, v, t) is computed at each time step on a fixed, uniform grid sampling the simulated phase space (x, v)). This approach has the merit of avoiding any statistical noise that is typical of particleincell (PIC) methods, but at the price of a large computational cost and memory requirements (when compared, for instance, to typical stateoftheart hybridPIC simulations with similar parameters; e.g., Cerri et al. 2017). To advance in time the proton distribution function, the HVM code implements an extended version of Knorr’s splitting scheme that accounts for the electromagnetic case (Mangeney et al. 2002). Furthermore, to efficiently couple the Vlasov equation for the proton distribution function with the Maxwell equations, the HVM code also employs the socalled current advancement method (CAM) (Matthews 1994). In the CAM method, in order to advance the distribution function from t to t + Δt, electromagnetic fields at the halfstep t + Δt/2 are needed. These fields depend upon the moments of the distribution function via the generalized Ohm’s law (3), so that in principle f_{p}(t + Δt/2) would also be needed. However, at this step of the algorithm, only f_{p}(t) and its “advectiononly advancement”^{1} are known from the splitting scheme. The CAM approach can use these known states to provide an approximation of the electromagnetic fields, which is sufficiently good in the sense that it does not lower the whole algorithm accuracy (which is secondorder in time). Further details about the HVM numerical scheme can be found in Valentini et al. (2007) and references therein.
2.1.1. HVM: Isothermal electron fluid
The simplest class of closures is represented by a isotropic, barotropic fluid, that is, where the pressure tensor is isotropic, Π_{e} = p_{e}I (I being the identity tensor), and the scalar pressure is a function of the density only, p_{e} = p_{e}(n_{e}). The first closure considered in our hybrid model is represented by isothermal electrons:
where T_{0, e} is the constant, homogeneous electron temperature.
2.1.2. HVLF: Landaufluid electrons
In magnetized, collisionless plasmas, the isotropy condition is rarely met. Charged particles indeed freely stream along B, while they gyrate around the field lines. It is therefore natural to distinguish between the parallel and perpendicular response of the plasma, especially for what concerns the distribution of energy between the different degrees of freedom. To the lowest (zeroth) order in an FLR expansion of the pressure tensor, this can be taken into account by a “gyrotropic” pressure (Chew et al. 1956),
where bb is the dyadic tensor of the magneticfield unit vector, b = B/B. Then, we explicitly solve the dynamic equations for the evolution of p_{e, ∥} and p_{e, ⊥}. These are obtained via integration over v of Eq. (1) (but for electrons) multiplied by the dyadic (v − u_{e})(v − u_{e}), under the assumption of gyrotropic pressure (see Hunana et al. 2019b, and references therein), and they read
where q_{∥, e} and q_{⊥, e} are the two scalars determining the parallel component of a gyrotropic heatflux tensor (i.e., the flux of parallel and perpendicular thermal energy along the magnetic lines, respectively). In the hybridVlasov code, these equations are used to advance in time p_{∥, e} and p_{⊥, e} employing a thirdorder RungeKutta method. Also, a substepping approach is implemented for these equations (meaning that, while the proton distribution function is advanced from t to t + Δt by means of a single time step of “length” Δt, the electrons pressure is simultaneously advanced with n_{sub} substeps, each of length Δt/n_{sub}). We remark here that this work is focused on the physics that the HVLF model can reproduce when compared to a fullkinetic or an HVM approach (in the context of standard 2D MR), whereas a detailed technical paper focusing on the numerical method is currently in preparation.
At this point, we still need a closure for q_{∥, e} and q_{⊥, e}. For instance, in the ChewGoldbergerLow (CGL, Chew et al. 1956) doubleadiabatic model they are set to zero (hence the name). Here, we adopt an LF type of closure for the electron heatflux coefficients, that is, a closure that provides an approximated linear response based on kinetic theory. This closure was presented in Sulem & Passot (2015), and its derivation consists in substituting the plasma response function by threepole and onepole Padé approximants in the formulas for the heat fluxes and the temperatures, respectively, which in turn are given by the linear kinetic theory. The closure is then expressed as closed formulas for the heatflux density coefficients, which in the case of the electrons are
Here T_{∥e} and T_{⊥e} denote the parallel and perpendicular electron temperatures, is the electron parallel thermal velocity of the unperturbed plasma (with m_{e} the electron mass), and Ω_{p} = eB_{0}/m_{p}c is the proton cyclotron frequency (with B_{0} the intensity of the magnetic field in the unperturbed state). In Eq. () the overbar, , denotes the space average and ℋ is the operator represented, in Fourier space, by
where the hat, , denotes Fouriertransformed quantities. This operator effectively takes into account the distortion of magneticfield lines in the nonlinear regime. This operator should have been the (negative) Hilbert operator along perturbed field lines going through each point, as in Eq. (50) of Snyder et al. (1997). This however involves nonaffordable computations in large simulations. Therefore, we adopted the semiphenomenological modeling of ℋ, still local in Fourier space, presented by Passot et al. (2014). In the linear regime, when b = e_{z}, the Fourier representation of ℋ given by Eq. (8) reduces to i sgn(k_{z}) = i k_{z}/k_{z}, as originally introduced by Hammett & Perkins (1990). This semiphenomenological operator is formally a good approximation only where the magnetic field is not strongly curved, a condition that could be challenged by a nonhomogeneous equilibrium (as in the case of this work). We also stress that the LF closure defined through Eq. (8) breaks down in the presence of magnetic nulls, since b cannot be defined where the magnetic field vanishes. In fact, the b ⋅ ∇ operator becomes illbehaved whenever B approaches zero. For this reason, in this work, we adopt a finite value for the initial guide field (low enough not to slow down the reconnection dynamics, but sufficiently large to avoid the emergence of magnetic nulls throughout the simulation).
Additionally, we mention that a possible approach could be to add also some artificial constraints to prevent runaway anomalies. This would be the case when hardcoded “anisotropy limiters” are employed to bound the values of the pressure anisotropy within a physically “sensible” range, something usually adopted when modeling (nearly) collisionless plasmas as fluids. For example, to mimic microphysical instabilities that “feed” on the pressure anisotropy, limiters were added to collisionless MHDlike models to perform numerical studies of magnetorotational instability (Sharma et al. 2006), to perform numerical studies of the turbulent dynamo in the intracluster medium (SantosLima et al. 2014; StOnge et al. 2020), and to derive a heating mechanism able to stably balance radiative cooling in the central regions of coldcore clusters of galaxies (Kunz et al. 2011). In all these examples the inclusion (via the limiters) of microphysical effects significantly changes the saturation level of macroscopic quantities. However, the current work seems to suggest that, in the present application, the implemented LF closure alone performs quite adequately to that end (as long as the conditions for the LF model to be valid are satisfied). From the figures presented in Sect. 3.2, one can see that the electron anisotropy remains always bounded within a reasonable range. Nevertheless, we do not exclude that including such limiters for a certain class of problems could be an interesting option to be further investigated in the future.
Finally, we note that in this model we do not include electrons’ FLR corrections to the gyrotropic pressure tensor in (5). Therefore, although retaining a model for the electron Landau damping, we do not capture other kinetic effects at k_{⊥}ρ_{e} ∼ 1, where ρ_{e} is the electron Larmor radius (see, e.g., Sulem & Passot 2015 for a detailed discussion). On the other hand, we stress that the electron model adopted here includes several additional effects with respect to the electron response adopted in a similar study by Le et al. (2016). In fact, while their model only considers a static, gyrotropic equation of state (viz., p_{∥, e} = p_{∥, e}(n, B) and p_{⊥, e} = p_{⊥, e}(n, B)), our HVLF model includes dynamic equations (6a)–(6b) for the gyrotropic electronpressure components, LF closures (7a)–(7b) for the gyrotropic electron heat fluxes, and finite electroninertia effect in the generalized Ohm’s law (3).
2.2. Fullkinetic model: iPIC
Compared to the LF closure implemented in the HVLF run, a fullkinetic model will provide also the nongyrotropic effects of the dynamics of the electrons (i.e., those linked with FLR physics). Also, the quasineutrality ansatz will be removed, meaning that charge separation is not forced to be zero.
The most widely used method to numerically describe the kinetic dynamics of a plasma system is the PIC method. The majority of PIC codes use an explicit method to advance the electric and magnetic fields in time (Bowers et al. 2008). The explicitness of the method introduces strong restrictions regarding the length and timescales to resolve in the numerical integration: (i) the plasma Debye length, the shortest length scale in the system that represents the characteristic screening scale length of a single charge by collective plasma motions, must be resolved by the spatial grid; (ii) the electron plasma frequency has to be resolved temporally; (iii) a Courant condition on the speed of light limits even further the grid size. Implicit and semiimplicit PIC methods have been developed to reduce such constraints. They are not subject to the aforementioned restrictions and allows for the simulation of fully kinetic plasma systems in a range that goes from fluid scale to a fraction of the electron skin depth.
Considering the size and the resolution of the system presented here, we opted for the use of the semiimplicit method implemented on the code “iPic3D” (Markidis & Lapenta 2010; Lapenta 2012). The semiimplicitness of the method arises from an approximated calculation of the source terms (currents and density) at the iterative step “n+1,” which are computed by a Taylor expansion in time of the same quantities at the step “n.” As explained in Markidis & Lapenta (2010), the semiimplicit method has a reduced computational cost as compared with the implicit method and, clearly, with every explicit method. As a drawback, the semiimplicit method requires the following stability condition to be satisfied: v_{th, e}Δt/Δx < 1, where v_{th, e} is the electron thermal speed, Δt is the time step, and Δx is the grid step. The latter is always less restrictive than the Courant condition on the speed of light c for the explicit methods, which reads cΔt/Δx < 1. The code iPic3D solves the full VlasovMaxwell set of equations for multiple species, in this case, protons and electrons. The physics resolved in a particular simulation depends on the system resolution and on the time step. Contrary to explicit methods where all relevant spatial and temporal scales are resolved, the semiimplicit method can resolve phenomena that have a spatial scale larger than the grid size and a frequency higher than the inverse of the time step. Thus, the presence or absence of phenomena such as quasineutrality violation depends on the numerical resolution. In this study, the implicit nature of the code is effective in removing any wave generated by charge separation, so that the quasineutrality is mostly preserved.
2.3. Simulation setup
For each of the plasma models mentioned above, we performed a numerical integration using the same initial condition on a 2D spatial domain of size L_{x} = 24πd_{p} and L_{y} = 12πd_{p}, discretized with N_{x} × N_{y} = 1024 × 512 grid points. The velocity space domain in the HVM and HVLF hybrid simulations is bounded by v≤6.4v_{th, p} in each velocity direction, where is the proton thermal velocity associated with the initial proton temperature T_{0, p}. The velocity grid has been discretized by 51^{3} grid points. The iPIC (fully kinetic) run employs 14000 particle per cell (PPC), 7000 for each species. The hybrid simulations also include finite electroninertia effects, and in all three models we adopted a reduced mass ratio of m_{p}/m_{e} = 100.
The initial setup consists of a “doubleHarris (DH) sheet” configuration (Harris 1962) to accommodate for periodic boundary conditions:
where L_{1} = 0.85d_{p} and L_{2} = 1.7d_{p} are the two shear widths, y_{1} = L_{1}L_{y}/[2(L_{1} + L_{2})] and y_{2} = L_{y} − (L_{2}L_{y})/[2(L_{1} + L_{2})] their positions, and B_{0} = 1 is a constant value used for normalization. Unlike the original Harris’ equilibrium, we also included a “weak” guide field in order to have a nonvanishing magnetic field to allow for the LF closure. In other words, if the magnetic field vanishes, equations for parallel and perpendicular quantities become unjustifiable. We set the guide field value to a quarter of the (asymptotic) reconnecting field B_{x}, so as not to slow down reconnection significantly (see, e.g., Ricci et al. 2004; Daughton & Karimabadi 2005; Shi et al. 2020).
The particle distribution function of the species a = p, e is initialized consistently with the magneticfield configuration in Eq. (9):
where the density profiles of the “background” and sheets populations, namely n_{b}, n_{1}, and n_{2}, respectively, are given by
and n_{0} = 1 is a constant value used for normalization (also, in our units the Boltzmann constant is κ_{B} = 1). Here m_{a} and T_{0, a} are the mass and the initial (uniform) temperature of the species, its initial thermal velocity and V_{1, a} = −2T_{0, a}/(B_{0}L_{1}) and V_{2, a} = 2T_{0, a}/(B_{0}L_{2}) the velocities related to the current sheets. If one considers a single current sheet, L_{2} → +∞, then Eqs. (9)–(11) form a consistent kinetic equilibrium (see Allanson et al. 2017). Due to the nonlinear nature of the VlasovMaxwell system of equations, the superposition of two Harris’ sheets is no longer a Vlasov equilibrium. Nevertheless, with our choice of parameters the two sheets are well separated in space: at the peak of the first current sheet (n_{1} ≃ n_{0}), the second one is negligible (n_{2} ∼ 10^{−10}n_{0}), and vice versa. Moreover, the choice of L_{2} = 2L_{1} makes the second current sheet quiescent over the entire simulation time (i.e., no reconnection nor other instabilities develop).
For our simulations, we considered^{2}β_{p} = 1, T_{0, e}/T_{0, p} = 1/4, and T_{0, ∥, a} = T_{0, ⊥, a}. The choice of adopting colder electrons follows from the absence of electronFLR effects in the HVLF model. At the same time, the hybrid models employed in this work do include electroninertia effects associated with d_{e}. Therefore, in order to resolve the electron inertial length d_{e} while keeping the electron Larmor radius below the grid resolution, we adopted β_{e} < 1. The initial equilibrium is perturbed by long wavelength, random phase magneticfield fluctuations with 1 ≤ kd_{p} ≤ 9, where (all simulations employ the same phases for the initial perturbations). By initializing a superposition of different modes, we do not impose a single X point to emerge at a given position. Furthermore, in order to mitigate the effect of PIC noise, the maximum amplitude ϵ ≡ max(δB) of the initial perturbations used in hybrid simulations, namely ϵ_{hyb} = 0.01, is increased by a factor of three in the iPIC case, ϵ_{pic} = 0.03. Other characterizing parameters for the simulations are Δx = Δy = 0.074d_{p} = 0.74d_{e}, (with 3 electrons substeps), , and for the iPIC simulation only, also v_{A}/c = 0.01, ω_{P,e}/Ω_{e} = 10, λ_{D}/d_{e} = 0.035. Within our choice of parameters, in the asymptotic region of the Harris field the electron gyromotion is characterized by Ω_{e} = 100Ω_{p} and ρ_{e} = 0.035d_{p} (we recall that m_{p}/m_{e} = 100).
Given Δx and Δt^{(iPIC)} defined above, the time and spatial (asymptotic) resolution of the iPIC simulation in terms of electrons’ quantities is thus and Δx = 2ρ_{e} (corresponding to a frequency and wavenumber resolution of ω/Ω_{e} ≲ 15.7 and k_{⊥}ρ_{e} ≲ 1.6, respectively). As a result, the electron gyration is not exactly resolved in space where the magnetic field is stronger (asymptotic region), but it is wellresolved in time everywhere, with at least ≈30 time steps in the worstcase scenario (if τ_{e} = 2π/Ω_{e} is the gyroperiod, then Δt^{(iPIC)} = τ_{e}/(10π)). Therefore, phenomena at the electron cyclotron frequency can be described within the whole simulation domain. On the other hand, a limitation exists on the range of spatial scales that are resolved: although finitek_{⊥}ρ_{e} effects can be captured by the iPIC simulation, to a certain extent, it is true that these effects do not extend far below the electronLarmorradius scale (at least in the asymptotic region, where the magnetic field is stronger; electrons’ gyration is better resolved within the Harris sheet, where the magnetic field strength drops by a factor of up to four with respect to the asymptotic region).
3. Numerical results
In this section, we present an analysis of MR simulations within the different physical models introduced in Sect. 2 (namely HVM, HVLF, and iPIC). Our analysis focuses on the following aspects: (i) the overall reconnection dynamics, (ii) the morphology of the EDR and outflows, and (iii) the pressureanisotropy production and regulation.
In particular, we highlight the capability of the HVLF simulation of capturing part of the physics^{3} missing in the HVM one and playing an important role in the dynamics, in agreement with the full PIC approach. Of course, the richer electron physics visible within the full PIC model is missing in the HVLF simulation. Nevertheless, this does not prevent the hybrid model adopted in the HVLF run from describing the reconnection dynamics quite satisfactorily.
3.1. Reconnection dynamics
Reconnection rate. We characterize the global reconnection dynamics by means of (i) the reconnected flux and the associated reconnection rate (in other words, the amount of magnetic flux advected through the X point) and (ii) the physical structures formed during the reconnection and their nonlinear evolution. In 2D, the reconnected flux F is obtained by calculating the difference between the values of the magnetic flux function Ψ evaluated at the X point and at the adjacent Opoint as follows (see Yeates & Hornig 2011):
The associated reconnection rate, R, is obtained as the time derivative of F:
Although it can be shown that R is equivalent to the more commonly used parallel electric field evaluated at the X point, E_{∥}(XPoint), we point out that in PIC simulations the Ψ function is less affected by numerical noise than the electric field. Moreover, it does not require timeaveraging to smooth out the oscillations typically observed as a consequence of a slight readjustment to the initial setup, rather than to the reconnection process itself. In addition, the evolution of the reconnected flux allows the linear and nonlinear phases of the reconnection process among the different runs to be synchronized (and, thus, to be properly compared).
In Fig. 1 we show the reconnected flux F (left panel) and the associated reconnection rate R (right panel) obtained from the different simulations. From these plots we distinguish three consecutive phases: the linear stage of the instability (up to ), a subsequent nonlinear “superexponential” phase (up to ), and a final quasistationary fully nonlinear regime. The first phase presents the usual exponential growth of the modes, the second one is characterized by some modes growing faster than exp(γ_{lin, k}t) (hence the name superexponential or quasiexplosive, as defined in Ottaviani & Porcelli 1993, and Aydemir 1992), the third by a quasistationary reconnection rate.
Fig. 1. Reconnected flux, F (left), and reconnection rate, R (right), for the main X point in the simulations. Within each panel, an inset shows a zoomin on a specific stage of these curves (thickened on the main plot), shifted in time, for the sake of comparison. 
A qualitatively good agreement between the different simulations during the linear phase is observed, as shown in the left panel of Fig. 1 and its inset. However, a detailed, quantitative comparison of the modes’ growth rates during the linear phase is limited by the uncertainties associated with finitePPC noise. Thus, we tested the convergence of the PIC simulation, in particular during the nonlinear stage, using different numbers of macroparticles, ranging from 800 to 7000 PPC per species (the higher number is the value used for the simulation presented here). This nonlinear phase of MR is the one of interest here. We focus on the latetime evolution when the process develops enough anisotropy to affect the dynamics, and thus the “robustness” of the electrons LF model can also be tested.
A qualitative agreement between the three simulations (see Fig. 1, right panel and its inset) persists during the “superexponential” phase, while they begin to differ during the quasistationary regime. This feature is the consequence of the different response times of the electron dynamics associated with the three different models. During the quasisteady stage, the (stationary) reconnection rate is higher for the HVM simulation (around a value ∼0.06), while the HVLF one, which includes electron anisotropy, shows a better agreement with the fully kinetic case (both reaching a quasisteady reconnection rate value around ∼0.04). The same trend for the quasistationary reconnection rates (i.e., higher with isotropic and isothermal electrons, lower in the fullkinetic case) was also reported in Le et al. (2016). Although the stationary reconnection rates are smaller than the “usual” value of about 0.1 (Cassak et al. 2017), it is worth noticing that this rate can indeed also depend on the background density (see Wu et al. 2011; Divin et al. 2019) or, more precisely, on the ratio between the background (i.e., asymptotic) and the peak (i.e., in the middle of the Harris sheet) values of the density. In fact, most simulations adopt a backgroundtopeak density ratio of n_{b}/n_{H, 0} = 0.1 − 0.2^{4}, in order to better compare with early simulations setup (usually GEMlike configurations, see Birn et al. 2001) and with some magnetospheric data. In Cassak et al. (2017) the socalled 0.1 reconnection rate problem is addressed showing that indeed not all in situ observations and/or simulations report exactly that value of the normalized reconnection rate. In our simulation, for instance, the backgroundtopeak density ratio is n_{b}/n_{1}(y_{1}) = 1.25 (cf. Eq. (11)) and thus a lower stationary reconnection rate is consistent with previous works (Divin et al. 2019).
Morphology. We now discuss the main features of the reconnection regions. In Figs. 2 and 3 we show for each simulation four fields of interest, namely ΔB_{z} = B_{z} − B_{z}(t = 0), E_{∥}, J_{z}, and , in the region around the X point during the nonlinear phase. The two figures are taken before and during the stationary phase, respectively. If not otherwise stated when comparing results from different simulations, the time is chosen such that the corresponding value of the reconnected flux F at that time is the same for all simulations (cf. Fig. 1). First of all, from Figs. 2 and 3 one observe that the morphology of the X point, as well as the symmetry of the various physical quantities around that region, highlights a significant difference between the three models, marking the distinction between the dynamics resulting from isothermal electrons (HVM) and from the electron’s models adopted in the HVLF and iPIC simulations. Moreover, as seen in the bottom row of these figures, in the HVLF and iPIC case one observes (more and more) elongated electron jets as compared to the HVM case. These jets, as well as the outofplane current, are tilted with respect to the direction of the neutral line. Another feature present in the HVLF and iPIC cases only is an asymmetric pattern in B_{z} fluctuations within the EDR. In the E_{∥} field, more clearly seen in Fig. 3, a similar structure is present in the three simulations. However, its typical extension is comparable in the HVLF and iPIC cases only, both quite larger than in the HVM one, during the quasistationary phase. In the PIC simulation, as compared with the other fields, E_{∥} is noisier (the reason why we prefer Ψ to evaluate F) and “fragmented.” This fragmentation, also visible in other quantities, is due to the formation ad absorption of a plasmoid in the time range from to (cf. Fig. 1). In this paper, we highlight other effects of this transient plasmoid; however, we can state that its formation did not pose a problem to the present study.
Fig. 2. Heat maps of ΔB_{z} = B_{z} − B_{z}(t = 0) (first row), E_{∥} (second row), J_{z} (third row), and (fourth row) for the HVM (left column, ), HVLF (central column, ), and iPIC (right column, ) simulations. These times are chosen by fixing the value of the reconnected magnetic flux. This value is in the middle of the superexponential phase. The plot is zoomedin around the main X point. The is the rms value of E_{∥} in the shown region. 
All the abovementioned features have already been pointed out by Egedal et al. (2009) as typical signatures for fully kinetic models in the presence of a guide field. The same authors proposed a fluid closure for electrons (see Le et al. 2009), which can be interpreted as an “interpolation” between an isothermal and a doubleadiabatic closure, with the addition of kinetic features^{5}. This closure has been validated in Le et al. (2009, 2010b) and Egedal et al. (2013), within a twofluid framework, and an analogous configuration for the outofplane electrons current was obtained. However, the current sheet in the simulation with their model forms a smaller angle with the neutral line, if compared with an analogous fullkinetic run. In this sense, the HVLF results are closer to those observed in the fullkinetic iPIC case.
Electronanisotropy effects on EDR shape. The morphology of the region around the X point is a consequence of the EDR elongation that arises both in the iPIC and HVLF simulations, but not in the HVM case. Although such an elongation is more pronounced in the iPIC simulation than it is for the HVLF one, this supports the idea that the elongation of the EDR is physical and associated with the local electron dynamics, the main reason being the electron pressure anisotropy. It is remarkable that a long time ago Vasyliunas (1975) and, twenty years later, Cai & Lee (1997), pointed out that, during steadystate MR, the forces associated with the anisotropy in the electron pressure should play a significant role in balancing the electric field associated with reconnection in the neighborhood of the X point. This idea was investigated more recently by Le et al. (2009); Le et al. (2010a,b) by means of a selfconsistent model involving electron pressure anisotropies, Hall magnetic field, and the electrons jets. Inside and near the EDR, as shown in Fig. 3, a parallel electric field is generated as a consequence of the breakdown of the frozenin condition. This electric field enables efficient parallel electron heating and, together with a simultaneous adiabatic perpendicular cooling (due to gradients in the magnetic field), produces a significant amount of electron anisotropy, p_{‖,e} ≫ p_{⊥,e}, typically observed in satellite data and confirmed in fully kinetic simulations. In our simulations, this anisotropy can be seen in the second row of Fig. 4. Here, we observe a difference between the iPIC and the HVLF simulations. In particular, in the iPIC simulation the electron anisotropy affects a larger area. This is consistent with the idea of heating driven by electron trapping (see Egedal et al. 2013, and references therein). The kinetic trapping process is obviously completely lost in the electron isothermal HVM model but also in the extended HVLF description because of the fluid approach. On the other hand, it is remarkable that the anisotropy level reached in the HVLF simulation is comparable to that measured in the iPIC simulation. This is in agreement with the observed larger E_{∥} structure in the HVLF and iPIC simulations. It is then the anisotropy of the electrons which in turn helps to balance the enhanced parallel electric field in the generalized Ohm’s Law.
Fig. 3. Same as Fig. 2, but the common value of reconnected flux is chosen in the nonlinear phase (for HVM at , for HVLF at , and for iPIC at ). 
Effects of electron anisotropy can be seen also in the momentum balance. Neglecting electron inertia and assuming a gyrotropic Π_{e}, electrons momentum balance requires a perpendicular current given by
The first term on the righthand side is the E × Bdrift. The two other drifts are driven by the pressuregradient term ∇ ⋅ Π_{e} (namely, the diamagnetic and the curvature drifts, respectively). In the region around the X point where electron jets develop, the last term of Eq. (14) dominates over the E × Bdrift due to the high anisotropy and strong curvature. In Fig. 4 the curvature b ⋅ ∇b is shown in the third row, while the ∇ ⋅ Π_{e}drift (J_{Π, e}) is shown in the fourth one. J_{Π, e} turns out to be much larger in the HVLF and iPIC simulations as compared with the HVM one, and since the curvature does not show significant differences between the three runs, the difference can be only due to the electron anisotropy. As shown in the figure, this J_{Π, e} driftcurrent has clearly the greatest influence on (cf. Fig. 3).
Fig. 4. Heat maps of proton anisotropy (first row), electron anisotropy (second row), magnetic curvature (third row), and ∇ ⋅ Π_{e}drift current (fourth row) for the HVM (left column, ), HVLF (center column, ), and iPIC (right column, ) simulations. The three instants are chosen by fixing the value of the reconnected flux. This value is chosen in the quasistationary phase. The plot is zoomedin around the main X point. 
This resulting anisotropydominated current is thus the main responsible for the asymmetric Hall magnetic field B_{H} observed in HVLF and iPIC simulations. Another consequence is that the J × B force becomes dominant^{6}. In order to understand what balances this force, the steadystate electron momentum equation can be cast in the form
with F representing the other contributions (electric field, inertia, nongyrotropy). Where the magneticfield curvature becomes important (also thanks to B_{H}), the term resulting from ∇ ⋅ (bb) can only be balanced by decreasing the quantity p_{∥, e} − p_{⊥, e} − B^{2} (i.e., by increasing p_{∥, e} − p_{⊥, e}). Thus the electron pressure anisotropy acts as to roughly balance the magnetic tension, which is the equivalent of saying that J × B ≈ ∇ ⋅ Π_{e}. To summarize, anisotropy and Bcurvature enhance electrons outflow jets, which generate the Hall magnetic field, B_{H}. This Hall component simultaneously provides the major contribution to the magneticfield curvature, which in turn requires strong anisotropy to achieve force balance. This seems to explain why extended narrow electrons jets forms in models that include the electrons pressure anisotropy.
3.2. Anisotropy
During the evolution of MR, we expect anisotropy to develop and to play a relevant role not only in the dynamics and morphology of the EDR around the X point (as outlined at the end of the previous subsection) but also by representing a source of free energy for possible secondary instabilities.
In Fig. 5 we show an histogram of the plasma distribution in a parameter space described by pressure anisotropy, A_{a} = p_{⊥, a}/p_{∥, a}, versus parallel plasma beta, β_{∥, a} = 8πp_{∥, a}/B^{2}, with a = p, e. Such a distribution is obtained by considering only the plasma within a region of the simulation that includes the active current sheet (roughly one third of the simulation box, in height). Also, we take the cumulative data over a time interval in the quasistationary phase^{7}. This is a typical representation of the plasma distribution used, for instance, by SW turbulence studies (e.g., Hellinger et al. 2006; Matteini et al. 2007; Bale et al. 2009). In these plots, we also report a set of curves representing the thresholds for anisotropydriven instabilities, listed in the top right box. The upper branches (p_{⊥, a} > p_{∥, a} region) are the MI (red) and the cyclotron instability (CI; orange); lower branches are parallel FHI (PFHI; green) and oblique FHI (OFHI; violet). The same color is used for both the protons and the electrons versions. These curves, used for reference, are taken from different works. For what concerns protons, the CI threshold is from Lazar & Poedts (2014), while the MI and FHI curves are taken from Maruca et al. (2012) and Astfalk & Jenko (2016), respectively. The electron version for these thresholds are instead taken from Lazar et al. (2013) (CI), Gary & Karimabadi (2006) (MI), Gary & Nishimura (2003) (PFHI), and Hellinger et al. (2014) (OFHI). In order to be consistent with the duration of our simulations, all the protoninstability thresholds, as well as the electron PFHI curve, are computed for the maximum growth rate . The remaining electroninstability thresholds are computed for (which, in our simulations, means ).
Fig. 5. Cumulative histogram of the (β_{∥, a}, T_{⊥, a}/T_{∥, a}) occurrences over the quasistationary phase, for protons (left) and electrons (right) in the HVM (top, ), HVLF (middle, ), and iPIC (bottom, ) simulations. The spatial region considered is defined by y ≤ 12.5d_{p}. Solid curves represent different instability thresholds (the same color is adopted for both the proton and electron versions): cyclotron instabilities (orange), mirror instabilities (red), parallel firehose instabilities (green), and oblique firehose instabilities (violet). The black dashdotted line is 1/β_{∥, a}. 
Proton anisotropy. By comparing the distribution of protons in HVM, HVLF, and iPIC simulations in Fig. 5, one can see that they do not behave qualitatively very differently. That is true also in the spatial domain, as presented in Fig. 4. In all the three cases (in Fig. 5), a significant fraction of protons are indeed exceeding the FHI thresholds. However, a difference in proton anisotropy can be seen in the HVM simulation, where there is a larger fraction of protons that reaches values of pressure anisotropy that are further below 1, with respect to the HVLF and iPIC cases. In fact, when electrons are isothermal, the free energy available in the system can only feed the protons’ anisotropy. On the other hand, when electron anisotropy is allowed, such free energy is shared between the two species. This suggests that the electron model implemented in the HVLF simulation is capable of distributing free energy among the two species in a way similar to the fullkinetic case.
The main question remains why a nonnegligible fraction of protons populate a region beyond the FHI marginal stability. In order to elucidate this point, we consider the spatial region confined within a magnetic island: analogously to Fig. 5, in Fig. 6 we show a timecumulative histogram where now only the plasma inside the magnetic island is considered (islands elements are selected using a threshold on the magnetic flux, namely Ψ/max(Ψ) > 0.9). The full circles reported on top of the plasma distribution represent a time series of the instantaneous values of anisotropy and parallel beta, averaged over the island domain, viz. (⟨p_{⊥, a}/p_{∥, a}⟩_{island}, ⟨β_{∥, a}⟩_{island}) versus time. The black circles denote the time range over which the cumulative histogram is computed, while green circles represent earlier simulation times. From Fig. 6, one clearly appreciates that the interior of the magnetic island (or, at least, part of it) is where the proton population that fills the FHIunstable regions comes from. On top of that, one clearly appreciates that there is a noticeably greater portion of protons that is beyond the FHI thresholds in the HVM simulation if compared with the HVLF and iPIC cases. In fact, there is a larger fraction of the proton distribution, especially in its lowβ_{∥, p} portion, that is found in the unstable region for the HVM case. This is highlighted also by the islandaveraged points, which, in the HVM simulation, reach the lowest values among the three models. These features of the proton anisotropybeta distribution within the various models is (one of) the effect of how adopting a different electron response feeds back onto the proton dynamics (namely, whether or not electrons can develop pressure anisotropy, and thus distribute the initial free energy among more available energy channels; see, e.g., Cerri & Camporeale 2020).
Fig. 6. Cumulative histogram of the (β_{∥, a}, T_{⊥, a}/T_{∥, a}) occurrences over the quasistationary phase for protons (left) and electrons (right) in the HVM (top, ), HVLF (middle, ), and iPIC (bottom, ) simulations. The spatial region considered is the island interior, here defined by Ψ/max(Ψ) > 0.9. The full circles on top of the histogram represent the time evolution of the islandaveraged anisotropy and parallel beta: black circles denote times over which the cumulative histogram is computed, while green circles represent earlier times. Curves are the same as in Fig. 5. 
Complementary to Fig. 6, we also performed an analysis of the spatial locations of the plasma elements out of thresholds within the simulation domain, using the whole distribution reported in Fig. 5 (not shown here). We find confirmation that the islands are the main origin of the proton population that is found in the parameter region unstable to the OFHI, although part of the outflows also contributes. At the same time, there are spatial regions that become unstable for the proton CI and which are similar in all three simulations: These regions are correlated with points where B_{y} accumulates due to the reconnection process. Finally, adjacent to these latter spatial domains, there are regions where protons are MIunstable.
The fact that the island is the unstable region is not surprising since two factors are simultaneously at play: (i) islands are located in a low magnetic field region, which favors the instability development (i.e., decreasing T_{⊥} due to magneticmoment conservation), and (ii) outflows and compression enhance parallel heating (i.e., increasing T_{∥}). The question is rather why this unstable region is not regulated by the instability (i.e., releasing anisotropyrelated freeenergy as, for instance, Alfvén waves, and bringing the distribution back to marginal stability). The FHI thresholds shown in the figures are computed for a (maximum) growth rate of 10^{−2}Ω_{p}, so that in about one should see, at least, a reduction of this unstable region. Our explanation is that the limited dimension of the islands (length < 40d_{p}) prevents the formation of modes at k ρ_{p} < 1, corresponding to the maximum growth rate for the FHI (Hellinger & Matsumoto 2000; Schekochihin et al. 2010). For and (see Fig. 6), one finds that k_{min} ρ_{p} ≈ 1. We conclude that the fastestgrowing modes associated with these FHI branches cannot emerge due to the spatial constraints imposed by the islands size.
In Fig. 6, dashdotted black lines are ∝1/β_{∥, a} curves, presented here to show how both electrons and protons populating the island seem to follow, on average, a doubleadiabatic evolution (i.e., close to a CGL model with zero heat fluxes; see Matteini et al. 2012; Hunana et al. 2019b), until . Of course, since proton FHI is geometrically suppressed, nothing confines them during this expansion, which continues until the quasistationary phase (i.e., the black circles). Nevertheless, in all three simulations, the protons average point seems unable to go back to marginal stability at later times.
Electron anisotropy. We now discuss the effectiveness of the HVLF model in capturing (as possible) the electrons dynamics in the (p_{⊥,e}/p_{‖,e},β_{‖,e}) plane, when compared to the fullkinetic case. In both HVLF and iPIC simulations, the electron distribution at p_{⊥}/p_{∥} < 1 is indeed well bounded by the FHI marginal stability thresholds (see Fig. 5). In particular, one can see that the distribution lower boundaries are slightly better in agreement with the OFHI than with the PFHI (besides, we remind the reader that the threshold for the parallel branch is computed for a maximum growth rate that is ten times smaller than its oblique counterpart). In Fig. 6, on the side of the electrons (right column), we can appreciate the confinement action of the electron OFHI within the island. Indeed, during its time evolution, the islandaveraged point roughly follows a doubleadiabatic trajectory, approaching the electrons FHI thresholds, and then it bounces back within the stable region. After that, from both the average points and the histograms, we see a tendency toward the electrons MI marginal stability threshold. A full (nonlinear) development of mirror modes should be associated with the formation of magnetic holes (or magnetic peaks, see Hellinger & Štverák 2018, and references therein). However, in our simulations, the presence of the MR process (and the associated magnetic structures) does not allow us to clearly identify the abovementioned mirrorrelated magnetic structures.
A similar degree of confinement is observed at high p_{⊥,e}/p_{‖,e} and high β_{‖,e} values, where the system is well restrained by the electrons MI for β_{‖,e} > 1. On the other hand, at lower β_{‖,e} values, where the electrons CI should be responsible for the confinement, the HVLF models (obviously) fails. This confinement is instead present in the fully kinetic model. Indeed, even if in Fig. 5 a small portion of the electrons are unstable also in the iPIC run, we clarify that this is due to a transient (namely, a plasmoid) observed in this simulation. This structure lasted for about , with a dynamics ten times faster than the one associated with the Xpoint reconnection. At this timescale (given the mass ratio implemented in our simulations), even the electrons “lag behind” the produced anisotropy. This anisotropy is quickly regulated after the plasmoid is “absorbed” within the main island. As a check for a link between this anisotropy and the plasmoid, we can show that it is absent in the island, as seen in Fig. 6. We also checked that the portion of electrons that are unstable to CI in the iPIC case is, indeed, spatially located within the plasmoid.
This result can be considered as a physiological difference between the electron models implemented in the HVLF and iPIC simulations (i.e., due to the lack of electronFLR effects in the HVLF model). On top of that, we note that the portion of electron distribution populating the CIunstable region in the HVLF simulation follows the 1/β_{∥, e} doubleadiabatic evolution (see Matteini et al. 2012; Hunana et al. 2019b). Moreover, we checked that this unstable portion of the electrons in the HVLF case is, again, inside the island. This can be viewed as one more signature that processes occurring within a magnetic island, such as island contraction, contribute significantly to the anisotropy generation in both directions (i.e., not only pushing the plasma toward FHI regimes). Overall, this corroborates the hypothesis that the differences in electrons CI confinement between HVLF and iPIC simulations are physical rather than numerical. We also mention that indeed the iPIC simulation presents a small amount of electron cooling due to the implicit algorithm: so, one may think that this cooling partially dissipates the freeenergy that would later lead to highp_{⊥,e}/p_{‖,e} anisotropies. Being aware of this numerical effect, we checked the behavior of internal, kinetic, and total energies (i) of the single species and (ii) of the whole plasma, during the iPIC simulation (not shown). Even if a certain amount of numerical cooling is occurring (mostly noticeable at early times, and at the expense of electrons internal energy), that amount is not enough to explain the differences just discussed between the iPIC and the HVLF runs, especially in the quasisteady stage. Therefore, we confirm that the lack of plasma population above the electrons CI threshold in the iPIC simulation is the result of a physical, rather than numerical, effect.
The iPIC and the HVLF models, as expected, support in some cases different dynamics, in particular for low β_{‖,e} (∼1) and high electron anisotropy (> 1). This discrepancy is due to the possibility of the iPIC model capturing electron anisotropylimiting processes that the HVLF cannot. Such kinds of processes can be reproduced only if the electron relevant scales (i.e., the electron Larmor radius in space and the electron gyration in time) are resolved in the model. In our iPIC simulation, we duly resolve electron gyration in time but not in space at the initial time (as discussed in Sect. 2.3). However, due to the system evolution, the electron gyroradius increases locally, becoming larger than the grid spacing in a few extended regions. Figure 7 shows a 2D map of the ratio between ρ_{e} and Δx, in which green and blue denote regions where this ratio is grater than one. Those regions present a higher perpendicular temperature that corresponds to a larger value of the local Larmor radius. It is in those regions that electron FLR effects are captured by the iPIC model. The consequence of these effects is to limit the further development of p_{⊥,e} > p_{‖,e} anisotropy, which, to the contrary, is not limited by these effects in the HVLF case. In conclusion, electron gyration scale effects are at play in the iPIC run and not in the HVLF one, resulting in different electron anisotropy dynamics in the two cases.
Fig. 7. Heat map of the local electron Larmor radius (normalized to the grid spacing) from the iPIC run, at . The color map is suited for a qualitative analysis: gray for ρ_{e}/Δx ≈ 1, red for ρ_{e}/Δx < 1, blue for 1 < ρ_{e}/Δx < 2, and green for ρ_{e}/Δx > 2. 
Character of magnetic fluctuations. In order to further investigate the role of these instabilities in the regulation of the electrons’ anisotropy distribution, we also analyze the properties of magnetic fluctuations in the quasistationary phase, in a similar way to what is done in Bale et al. (2009) for SW turbulent fluctuations. Such an analysis is reported in Fig. 8, where the properties of magneticfield fluctuations are reported in a (β_{∥, e}, p_{⊥, e}/p_{∥, e}) plane: this is done relating the fluctuations’ values at any spatial point to the anisotropy and parallel beta values of the plasma at that same position. (The curves and time intervals considered here are the same for Fig. 5, although a slightly more coarsegrained binning is employed for these plots.) The colormap refers to the binaveraged value of (top row), (middle row), or δB_{∥}/δB_{⊥} (bottom row), measured in the HVLF (left column) and iPIC (right column) simulations. Here, δB = δB≡B − ⟨B⟩_{t} and , where ⟨…⟩_{t} represents time average over the range that is being considered. If just the magnitude of magneticfield fluctuations are considered (i.e., ; Fig. 8, top row), one finds that the larger fluctuations are mainly located in the regions nearby marginalstability thresholds, for both p_{⊥, e}/p_{∥, e} > 1 and p_{⊥, e}/p_{∥, e} < 1 values of anisotropy. Such a fluctuation enhancement is interpreted as a signature of the growth of these anisotropydriven instabilities in their attempt to regulate electronpressure anisotropy: larger fluctuations tend to scatter electrons more efficiently, so to isotropize their pressure and bring the plasma back to marginal stability. Considering now the magnetic fluctuations transverse to the mean field, (Fig. 8, middle row), one finds that these type of fluctuations are strongly enhanced across the FHI marginal stability (at p_{⊥, e}/p_{∥, e} < 1), and nearby the MI and CI thresholds (at p_{⊥, e}/p_{∥, e} < 1), but only for those regions where β_{∥, e} ≳ 10). These features are slightly clearer in the iPIC simulation than they are in the HVLF case. Indeed, in the HVLF simulation, the electrons do not clearly reach and/or overcome the FHI thresholds at (β_{∥, e} < 10). As mentioned above, this difference is due to the presence of a transient plasmoid structure in the iPIC simulation (cf. Figs. 5 and 6).
Fig. 8. Cumulative, binaveraged histogram of (top row), (middle row), and δB_{∥}/δB_{⊥} (bottom row) in the (β_{∥, e}, T_{⊥, e}/T_{∥, e}) plane. Histograms are cumulated over the quasistationary phase ( for HVLF, left column; for iPIC, right column). Here, δB = δB≡B − ⟨B⟩_{t} and , where ⟨…⟩_{t} denotes the time average over the interval used for the histogram. The spatial region considered and the curves are the same as in Fig. 5. 
Considering now the socalled magnetic compressibility of the fluctuations, δB_{∥}/δB_{⊥} (Fig. 8, bottom row), both simulations exhibit an enhancements of this fluctuations’ property near the MI (and/or CI) thresholds, when 1 ≲ β_{∥, e} ≲ 10. Magneticfluctuations compressibility is also found to be enhanced near the FHI thresholds, and in particular around the oblique branch, for those regions where 5 ≲ β_{∥, e} ≲ 100. What is often called OFHI should be better called nonpropagating FHI, since at oblique propagation the firehose splits into a propagating, ω_{r} ≠ 0, and a nonpropagating, ω_{r} = 0, branches (here, ω_{r} is the realpart of the complex frequency of the mode, ω = ω_{r} + iγ). The ω_{r} ≠ 0 branch is a continuation (at nonzero propagation angles) of the parallel case, it is a transverse (noncompressional) mode, and it resonates with protons, not with electrons. The nonpropagating branch instead is (at least partially) compressive and dominates over the propagating one. Focusing on the ω_{r} = 0 branch, when the compressive component is small (i.e., δB_{∥}/δB_{⊥} < 1), the instability is cyclotron resonant with electrons, and for this reason is able to control the electron anisotropy. When δB_{∥}/δB_{⊥} > 1, this mode becomes an analogous to the electronmirror mode, just on the other “side” of the anisotropy spectrum (see Lazar et al. 2014, and references therein).
Of course, it is important to stress also the limits of the HVLF model. We already mentioned the inability of the LF closure to reproduce the electrons CI. In the bottom row of Fig. 8, from the iPIC simulation we observe a very small magnetic compressibility for β_{∥, e} < 1, especially near the instability thresholds at p_{⊥, e}/p_{∥, e} > 1. This is likely due to the fact that electrons anisotropy is efficiently confined by their CI in the iPIC simulation. On the other hand, it is clear that the HVLF model cannot account for the electron CI process: the high magneticcompressibility values in this region thus suggest that is the MI that tries to make up for the absence of the CI in regulating the electronpressure anisotropy (with little success).
4. Discussion and conclusions
In this study we performed 2D numerical simulations of MR within a typical Harris’s sheet configuration. The same initial configuration is evolved by means of three different plasma models: (i) an HVM model with (isotropic) isothermal electrons with mass, (ii) an HVLF model, where an anisotropic electron fluid is equipped with an LF closure, and (iii) a fullkinetic model employing an iPIC algorithm. The goal of this work is to investigate to which extent different reduced electron models are able to capture the physics of MR, with the focus on its nonlinear stage. In particular, we are interested in elucidating the advances represented by the HVLF model, with respect to a standard HVM approach, in the capability of capturing several aspects of the electron dynamics and its feedback onto the protons. These aspects include the development of electronpressure anisotropy and anisotropydriven (electron) instabilities, as well as the electronkinetic response, such as Landau damping. Limitations that still hold within an HVLF approach with respect to a fullkinetic model are also highlighted.
As far as the global reconnection dynamics, expressed by the reconnected flux and its reconnection rate, is concerned, all three simulations show results that are qualitatively in agreement with one another. In particular, after an initial linear stage exhibiting minor differences between the HVM model and the other two models, the three simulations achieve a very good agreement in the socalled superexponential phase. In the quasistationary stage, on the other hand, there is a quantitative difference in the quasisteady reconnection rate between HVM and the other two models that is worth mentioning: The reconnection rate in the quasisteady regime is noticeably higher in the HVM model, about ≈0.06, than in the HVLF and iPIC simulations, where it sets around ≈0.04. This reconnection rate value differs from the one considered to be a “typical” value for fast collisionless reconnection, ≈0.1. Such a difference in the quasisteady reconnection rate is due to a backgroundtopeak density ratio that in our configuration is higher than the one that is “typically” implemented (see, e.g., Cassak et al. 2017, and references therein). As for the local features of the dynamics, the absence of electron anisotropy in the HVM simulation results in a different shape of the region around the X point of the reconnecting current sheet and, in particular, of the EDR. From the X point, in a region surrounding the EDR, elongated and “tilted” electron jets develop as a consequence of the electronpressure anisotropy that affects the momentum balance in the HVLF and iPIC simulations. For the same reason, the current structure and the EDR are also elongated (and the former is also tilted), a feature that is not observed with isotropic isothermal electrons.
Then, we focused on the production and regulation of species’ pressure anisotropy during the nonlinear stage of reconnection. The processes mainly involved in the confinement of the plasma distribution in a anisotropyversusbeta parameter space are the OFHI (for p_{⊥, a} < p_{∥, a}), the MI (for p_{⊥, a} > p_{∥, a} and β_{∥, a} ≳ 2), and the CI (for p_{⊥, a} > p_{∥, a} and β_{∥, a} ≲ 2). As far as the dynamics of the protons is concerned, all three simulations are qualitatively similar. A significant portion of protons that are inside a magnetic island continuously populate the parameterspace region that is beyond the OFHI marginalstability threshold. This is due to the size of the island not being long enough for the fastest growing mode of the instability to develop and to efficiently regulate their anisotropy (while island contraction simultaneously keeps generating more pressure anisotropy). Furthermore, when compared to the HVLF and iPIC simulations, protons in the HVM simulation show a tendency to develop larger values of pressure anisotropy (toward p_{⊥, i} < p_{∥, i}). This is the consequence in the HVM model of adopting an isotropic pressure for the electrons, which excludes a priori the possibility of feeding electron anisotropy using part of the energy released by the reconnection process. However, we defer a detailed study of this aspect to future works.
As far as electron dynamics is concerned, we have shown that the HVLF and iPIC simulations lead to very similar results, although some differences between the two are still present. A key finding is that, when compared to the fullkinetic case, the anisotropic electronfluid model with an LF closure efficiently regulates the distribution of electrons in an anisotropyversusbeta plane. In fact, the OFHI and the MI are effective at controlling the electron anisotropy in both simulations. The main difference between HVLF and iPIC arises due to the intrinsic limitation of the electronLF model, which cannot capture the electron CI. The CI threshold is indeed effective at limiting the electronpressure anisotropy (at p_{⊥, e} > p_{∥, e} and β_{∥, e} ≲ 2) in the iPIC simulation, while in the HVLF case such anisotropy can only be regulated by the MI (which in our setting is less effective than CI, and thus a larger electron anisotropy is achieved). This means that, when electrongyration effects are taken into account, such large anisotropy is not observed. In this context, electron FLR effects could also play a relevant role in limiting the anisotropy.
In conclusion, we have employed three different models to describe the same process, focusing on the effectiveness of how the electron physics can be modeled within a hybridkinetic framework (namely, including electronpressure anisotropy and a fluid model of their linear kinetic response). When compared to the fullkinetic case, the electron response implemented in the HVLF model effectively reproduces the main features of MR in 2D, as well as several aspects of the associated electron microphysics and its feedback onto protons dynamics. This includes the global evolution of MR and the local physics occurring within the EDR, as well as the evolution of species’ pressure anisotropy. In particular, OFHI, MI, and CI play a relevant role in regulating electrons’ anisotropy during the nonlinear stage of MR. As expected, the HVLF model captures all these features, except for the electron CI. The HVLF model represents a good approximation of the fullkinetic case that compromises between the richness of the physics that is described and its computational cost.
Future works will extend this study to MR in 3D, as well as to plasma turbulence across the socalled transition range between MHD and subion scales (where reconnection plays a crucial role; see, e.g., Cerri & Califano 2017, and references therein). Regarding additional model development, one can further refine the electronfluid model by including (largescale) electron FLR effects and/or moving the LF closure to higher moments (see, e.g., Sulem & Passot 2015). This would make the electron model even closer to a fullkinetic response, although at a somewhat higher computational cost, and is currently considered to be a possible future direction.
Actually, it is derived from considerations about the distribution function, discriminating between trapped and passing electrons. Thus, it is more than a naive interpolation between two preexisting models. However, in the weaktrapping and strongtrapping limits it approaches the isothermal and doubleadiabatic models, respectively.
Equation (14) implies . When comparing the Hall force associated with such a current, , with the force associated with the reconnecting electric field, n e E_{rec} ∼ n e (0.1 v_{A} B), one finds that J × B≫n e E_{rec} for the typical values of β_{e} investigated in the present work.
Since we “synchronize” the simulations via the reconnected magnetic flux, these intervals differ from run to run and are reported in the description of Fig. 5.
Acknowledgments
This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 776262 (AIDA, www.aidaspace.eu). F. F., S. S. C., and F. C. want to thank D. Telloni for useful discussions regarding insitu observations associated with magnetic reconnection. We acknowledge the Italian ISCRA initiative (grant HP10B10ALD), the European PRACE initiative (grant n. 2017174107 and n. pn56ye) for awarding us access to the supercomputer Marconi, CINECA, Italy and to SuperMUCNG at GCS@LRZ, Germany. We thank Dr. M. Guarrasi (CINECA, Italy) for his essential contribution to code implementation on Marconi. S.S.C. is supported by the MaxPlanck/Princeton Center for Plasma Physics (NSF grant PHY1804048). F.P. is supported by the PostDoctoral Fellowship 12X0319N and the Research Grant 1507820N from Research Foundation – Flanders (FWO). This work has been performed under the Project HPCEUROPA3 (INFRAIA20161730897), with the support of the EC Research Innovation Action under the H2020 Programme; in particular, F. P. acknowledges the support of F. C. and The Department of Physics at the University of Pisa and the computer resources and technical support provided by CINECA.
References
 Allanson, O., Wilson, F., Neukirch, T., Liu, Y. H., & Hodgson, J. D. B. 2017, Geophys. Res. Lett., 44, 8685 [NASA ADS] [CrossRef] [Google Scholar]
 Astfalk, P., & Jenko, F. 2016, J. Geophys. Res. (Space Phys.), 121, 2842 [NASA ADS] [CrossRef] [Google Scholar]
 Aydemir, A. Y. 1992, Physics of Fluids B, 4, 3469 [NASA ADS] [CrossRef] [Google Scholar]
 Bale, S. D., Kasper, J. C., Howes, G. G., et al. 2009, Phys. Rev. Lett., 103, 211101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Birn, J., Drake, J. F., Shay, M. A., et al. 2001, J. Geophys. Res., 106, 3715 [Google Scholar]
 Biskamp, D. 1994, Phys. Rep., 237, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Bowers, K. J., Albright, B., Yin, L., Bergen, B., & Kwan, T. 2008, Phys. Plasmas, 15, 055703 [NASA ADS] [CrossRef] [Google Scholar]
 Byers, J. A., Cohen, B. I., Condit, W. C., & Hanson, J. D. 1978, J. Comput. Phys., 27, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Cai, H. J., & Lee, L. C. 1997, Phys. Plasmas, 4, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Califano, F., Cerri, S. S., Faganello, M., et al. 2020, Front. Phys., 8, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Califano, F., Hellinger, P., Kuznetsov, E., et al. 2008, J. Geophys. Res. (Space Phys.), 113, A08219 [NASA ADS] [Google Scholar]
 Camporeale, E., & Burgess, D. 2017, J. Plasma Phys., 83, 535830201 [CrossRef] [Google Scholar]
 Caprioli, D., & Spitkovsky, A. 2013, ApJ, 765, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Cassak, P. A., Liu, Y. H., & Shay, M. A. 2017, J. Plasma Phys., 83, 715830501 [CrossRef] [Google Scholar]
 Cerri, S. S., & Califano, F. 2017, New J. Phys., 19, 025007 [NASA ADS] [CrossRef] [Google Scholar]
 Cerri, S. S., & Camporeale, E. 2020, Phys. Plasmas, 27, 082102 [CrossRef] [Google Scholar]
 Cerri, S. S., Franci, L., Califano, F., Landi, S., & Hellinger, P. 2017, J. Plasma Phys., 83, 705830202 [CrossRef] [Google Scholar]
 Cerri, S. S., Kunz, M. W., & Califano, F. 2018, ApJ, 856, L13 [CrossRef] [Google Scholar]
 Cerri, S. S., Grošelj, D., & Franci, L. 2019, Front. Astron. Space Sci., 6, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, L. J., Bessho, N., Lefebvre, B., et al. 2008, J. Geophys. Res. (Space Phys.), 113, A12213 [NASA ADS] [CrossRef] [Google Scholar]
 Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, Proc. R. Soc. London Ser. A, 236, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Comisso, L., & Sironi, L. 2018, Phys. Rev. Lett., 121, 255101 [NASA ADS] [CrossRef] [Google Scholar]
 Daughton, W., & Karimabadi, H. 2005, J. Geophys. Res. (Space Phys.), 110, A03217 [NASA ADS] [CrossRef] [Google Scholar]
 Daughton, W., & Karimabadi, H. 2007, Phys. Plasmas, 14, 072303 [Google Scholar]
 Divin, A., Semenov, V., Zaitsev, I., et al. 2019, Phys. Plasmas, 26, 102305 [NASA ADS] [CrossRef] [Google Scholar]
 Egedal, J., Le, A., & Daughton, W. 2013, Phys. Plasmas, 20, 061201 [NASA ADS] [CrossRef] [Google Scholar]
 Egedal, J., Daughton, W., Drake, J. F., Katz, N., & Lê, A. 2009, Phys. Plasmas, 16, 050701 [NASA ADS] [CrossRef] [Google Scholar]
 Franci, L., Cerri, S. S., Califano, F., et al. 2017, ApJ, 850, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Gary, S. P., & Karimabadi, H. 2006, J. Geophys. Res. (Space Phys.), 111, A11224 [NASA ADS] [CrossRef] [Google Scholar]
 Gary, S. P., & Nishimura, K. 2003, Phys. Plasmas, 10, 3571 [NASA ADS] [CrossRef] [Google Scholar]
 Grošelj, D., Cerri, S. S., Bañón Navarro, A., et al. 2017, ApJ, 847, 28 [CrossRef] [Google Scholar]
 Hammett, G. W., & Perkins, F. W. 1990, Phys. Rev. Lett., 64, 3019 [CrossRef] [PubMed] [Google Scholar]
 Harris, E. G. 1962, Il Nuovo Cimento, 23, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Hellinger, P., & Matsumoto, H. 2000, J. Geophys. Res., 105, 10519 [NASA ADS] [CrossRef] [Google Scholar]
 Hellinger, P., & Štverák, Š. 2018, J. Plasma Phys., 84, 905840402 [NASA ADS] [CrossRef] [Google Scholar]
 Hellinger, P., Trávníček, P., Kasper, J. C., & Lazarus, A. J. 2006, Geophys. Res. Lett., 33, L09101 [NASA ADS] [CrossRef] [Google Scholar]
 Hellinger, P., Trávníček, P. M., Decyk, V. K., & Schriver, D. 2014, J. Geophys. Res. (Space Phys.), 119, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Hellinger, P., Matteini, L., Landi, S., et al. 2019, ApJ, 883, 178 [NASA ADS] [CrossRef] [Google Scholar]
 Henri, P., Cerri, S. S., Califano, F., et al. 2013, Phys. Plasmas, 20, 102118a [NASA ADS] [CrossRef] [Google Scholar]
 Hewett, D. W., & Nielson, C. W. 1978, J. Comput. Phys., 29, 219 [CrossRef] [Google Scholar]
 Hunana, P., Tenerani, A., Zank, G. P., et al. 2019a, J. Plasma Phys., 85, 205850602 [CrossRef] [Google Scholar]
 Hunana, P., Tenerani, A., Zank, G. P., et al. 2019b, J. Plasma Phys., 85, 205850603 [CrossRef] [Google Scholar]
 Kunz, M. W., Schekochihin, A. A., Cowley, S. C., Binney, J. J., & Sanders, J. S. 2011, MNRAS, 410, 2446 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014, Phys. Rev. Lett., 112, 205003 [NASA ADS] [CrossRef] [Google Scholar]
 Lapenta, G. 2012, J. Comput. Phys., 231, 795 [CrossRef] [Google Scholar]
 Lazar, M., & Poedts, S. 2014, MNRAS, 437, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Lazar, M., Poedts, S., & Michno, M. J. 2013, A&A, 554, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lazar, M., Poedts, S., Schlickeiser, R., & Ibscher, D. 2014, Sol. Phys., 289, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Le, A., Egedal, J., Daughton, W., Fox, W., & Katz, N. 2009, Phys. Rev. Lett., 102, 085001 [NASA ADS] [CrossRef] [Google Scholar]
 Le, A., Daughton, W., Karimabadi, H., & Egedal, J. 2016, Phys. Plasmas, 23, 032114 [NASA ADS] [CrossRef] [Google Scholar]
 Le, A., Egedal, J., Fox, W., et al. 2010a, Phys. Plasmas, 17, 055703 [NASA ADS] [CrossRef] [Google Scholar]
 Le, A., Egedal, J., Daughton, W., et al. 2010b, Geophys. Res. Lett., 37, L03106 [NASA ADS] [CrossRef] [Google Scholar]
 Lembège, B., Savoini, P., Hellinger, P., & Trávníček, P. M. 2009, J. Geophys. Res. (Space Phys.), 114, A03217 [Google Scholar]
 Mangeney, A., Califano, F., Cavazzoni, C., & Travnicek, P. 2002, JCoPh, 179, 495 [NASA ADS] [Google Scholar]
 Markidis, S., Lapenta, G., et al. 2010, Math. Comput. Simul., 80, 1509 [CrossRef] [Google Scholar]
 Maruca, B. A., Kasper, J. C., & Gary, S. P. 2012, ApJ, 748, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Matteini, L., Landi, S., Hellinger, P., & Velli, M. 2006, J. Geophys. Res. (Space Phys.), 111, A10101 [NASA ADS] [CrossRef] [Google Scholar]
 Matteini, L., Landi, S., Hellinger, P., et al. 2007, Geophys. Res. Lett., 34, L20105 [Google Scholar]
 Matteini, L., Hellinger, P., Landi, S., Trávníček, P. M., & Velli, M. 2012, Space Sci. Rev., 172, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Matthews, A. P. 1994, J. Comput. Phys., 112, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, P. A., & Büchner, J. 2018, Phys. Rev. E, 98, 043205 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, P. A., Jain, N., Kilian, P., & Büchner, J. 2018, Comput. Phys. Commun., 224, 245 [CrossRef] [Google Scholar]
 Ottaviani, M., & Porcelli, F. 1993, Phys. Rev. Lett., 71, 3802 [NASA ADS] [CrossRef] [Google Scholar]
 Palmroth, M., Hoilijoki, S., Juusola, L., et al. 2017, Ann. Geophys., 35, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Passot, T., & Sulem, P. L. 2007, Phys. Plasmas, 14, 082502 [CrossRef] [Google Scholar]
 Passot, T., Henri, P., Laveder, D., & Sulem, P.L. 2014, EPJD, 68, 207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pucci, F., Servidio, S., SorrisoValvo, L., et al. 2017, ApJ, 841, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Ricci, P., Brackbill, J., Daughton, W., & Lapenta, G. 2004, Phys. Plasmas, 11, 4102 [NASA ADS] [CrossRef] [Google Scholar]
 Rincon, F., Califano, F., Schekochihin, A. A., & Valentini, F. 2016, Proc. Natl. Acad. Sci., 113, 3950 [NASA ADS] [CrossRef] [Google Scholar]
 SantosLima, R., de Gouveia Dal Pino, E. M., Kowal, G., et al. 2014, ApJ, 781, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Cowley, S. C., Dorland, W., et al. 2008, Plasma Phys. Controlled Fusion, 50, 124024 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Cowley, S. C., Rincon, F., & Rosin, M. S. 2010, MNRAS, 405, 291 [NASA ADS] [Google Scholar]
 Servidio, S., Dmitruk, P., Greco, A., et al. 2011, Nonlinear Processes Geophys., 18, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Servidio, S., Valentini, F., Perrone, D., et al. 2015, J. Plasma Phys., 81, 325810107 [CrossRef] [Google Scholar]
 Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. M. 2006, ApJ, 637, 952 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, C., Velli, M., Pucci, F., Tenerani, A., & Innocenti, M. E. 2020, ApJ, 902, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Snyder, P. B., Hammett, G. W., & Dorland, W. 1997, Phys. Plasmas, 4, 3974 [CrossRef] [Google Scholar]
 StOnge, D. A., & Kunz, M. W. 2018, ApJ, 863, L25 [CrossRef] [Google Scholar]
 StOnge, D. A., Kunz, M. W., Squire, J., & Schekochihin, A. A. 2020, J. Plasma Phys., 86, 905860503 [NASA ADS] [CrossRef] [Google Scholar]
 Sulem, P. L., & Passot, T. 2015, J. Plasma Phys., 81, 325810103 [CrossRef] [Google Scholar]
 Told, D., Cookmeyer, J., Muller, F., Astfalk, P., & Jenko, F. 2016, New J. Phys., 18, 065011 [NASA ADS] [CrossRef] [Google Scholar]
 Treumann, R. A., & Baumjohann, W. 2013, Front. Phys., 1, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Tronci, C., & Camporeale, E. 2015, Phys. Plasmas, 22, 020704 [NASA ADS] [CrossRef] [Google Scholar]
 Valentini, F., Trávníček, P., Califano, F., Hellinger, P., & Mangeney, A. 2007, JCoPh, 225, 753 [NASA ADS] [Google Scholar]
 Vasyliunas, V. M. 1975, Rev. Geophys. Space Phys., 13, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, H., Lin, Y., Wang, X., & Guo, Z. 2019, Phys. Plasmas, 26, 072102 [NASA ADS] [CrossRef] [Google Scholar]
 Weidl, M. S., Winske, D., Jenko, F., & Niemann, C. 2016, Phys. Plasmas, 23, 122102 [CrossRef] [Google Scholar]
 Winske, D. 1985, Space Sci. Rev., 42, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, P., Shay, M. A., Phan, T. D., Oieroset, M., & Oka, M. 2011, Phys. Plasmas, 18, 111204 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., & Hornig, G. 2011, Phys. Plasmas, 18, 102118 [CrossRef] [Google Scholar]
 Zhdankin, V., Uzdensky, D. A., Perez, J. C., & Boldyrev, S. 2013, ApJ, 771, 124 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1. Reconnected flux, F (left), and reconnection rate, R (right), for the main X point in the simulations. Within each panel, an inset shows a zoomin on a specific stage of these curves (thickened on the main plot), shifted in time, for the sake of comparison. 

In the text 
Fig. 2. Heat maps of ΔB_{z} = B_{z} − B_{z}(t = 0) (first row), E_{∥} (second row), J_{z} (third row), and (fourth row) for the HVM (left column, ), HVLF (central column, ), and iPIC (right column, ) simulations. These times are chosen by fixing the value of the reconnected magnetic flux. This value is in the middle of the superexponential phase. The plot is zoomedin around the main X point. The is the rms value of E_{∥} in the shown region. 

In the text 
Fig. 3. Same as Fig. 2, but the common value of reconnected flux is chosen in the nonlinear phase (for HVM at , for HVLF at , and for iPIC at ). 

In the text 
Fig. 4. Heat maps of proton anisotropy (first row), electron anisotropy (second row), magnetic curvature (third row), and ∇ ⋅ Π_{e}drift current (fourth row) for the HVM (left column, ), HVLF (center column, ), and iPIC (right column, ) simulations. The three instants are chosen by fixing the value of the reconnected flux. This value is chosen in the quasistationary phase. The plot is zoomedin around the main X point. 

In the text 
Fig. 5. Cumulative histogram of the (β_{∥, a}, T_{⊥, a}/T_{∥, a}) occurrences over the quasistationary phase, for protons (left) and electrons (right) in the HVM (top, ), HVLF (middle, ), and iPIC (bottom, ) simulations. The spatial region considered is defined by y ≤ 12.5d_{p}. Solid curves represent different instability thresholds (the same color is adopted for both the proton and electron versions): cyclotron instabilities (orange), mirror instabilities (red), parallel firehose instabilities (green), and oblique firehose instabilities (violet). The black dashdotted line is 1/β_{∥, a}. 

In the text 
Fig. 6. Cumulative histogram of the (β_{∥, a}, T_{⊥, a}/T_{∥, a}) occurrences over the quasistationary phase for protons (left) and electrons (right) in the HVM (top, ), HVLF (middle, ), and iPIC (bottom, ) simulations. The spatial region considered is the island interior, here defined by Ψ/max(Ψ) > 0.9. The full circles on top of the histogram represent the time evolution of the islandaveraged anisotropy and parallel beta: black circles denote times over which the cumulative histogram is computed, while green circles represent earlier times. Curves are the same as in Fig. 5. 

In the text 
Fig. 7. Heat map of the local electron Larmor radius (normalized to the grid spacing) from the iPIC run, at . The color map is suited for a qualitative analysis: gray for ρ_{e}/Δx ≈ 1, red for ρ_{e}/Δx < 1, blue for 1 < ρ_{e}/Δx < 2, and green for ρ_{e}/Δx > 2. 

In the text 
Fig. 8. Cumulative, binaveraged histogram of (top row), (middle row), and δB_{∥}/δB_{⊥} (bottom row) in the (β_{∥, e}, T_{⊥, e}/T_{∥, e}) plane. Histograms are cumulated over the quasistationary phase ( for HVLF, left column; for iPIC, right column). Here, δB = δB≡B − ⟨B⟩_{t} and , where ⟨…⟩_{t} denotes the time average over the interval used for the histogram. The spatial region considered and the curves are the same as in Fig. 5. 

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