© 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., fluid-scale) 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 magnetic-field topology of a system by locally breaking the magnetic connections ideally preserved at magnetohydrodynamic (MHD) scales (Biskamp 1994). The change in magnetic-field 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) “multi-scale 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, multi-scale 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 full-kinetic 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, full-kinetic 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 ideal-MHD theory represents the simplest case of an infinitely conducting, single-fluid 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 quasi-neutral 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 mass-less 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 mass-less 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 electron-pressure anisotropy (Chew et al. 1956; Hunana et al. 2019b), electron finite-Larmor-radius (FLR) effects, and/or electron Landau damping (e.g., Hammett & Perkins 1990; Sulem & Passot 2015; Hunana et al. 2019a).

In general, hybrid-kinetic 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; St-Onge & 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 kinetic-scale 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 hybrid-kinetic 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 magnetic-field direction, b = B/|B|), and a Landau-fluid (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–Landau-fluid” (HVLF). The idea is to include within a hybrid description the relevant electron pressure-anisotropy effects and a fluid model for the electron-kinetic 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 electron-kinetic effects at sub-ion 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 so-called 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 free-energy source for secondary instabilities, most notably fire-hose 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 electron-inertia effects and (ii) a fully kinetic model employing the semi-implicit particle-in-cell “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 hybrid-kinetic 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 full-kinetic model. It should be noted that the hybrid-kinetic approximation also keeps finite-electron-mass effects in both cases.

It is important to point out that, while classic hybrid-kinetic 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. Hybrid-kinetic models

The hybrid-kinetic approximations consists of a quasi-neutral plasma, n_p = n_e ≡ n, where fully kinetic ions (protons in this case) are coupled to a fluid-electron model (Winske 1985). The evolution of the proton distribution function, f_p(x, v, t), follows the Vlasov equation,

(1)

while the magnetic field, B, evolves following Faraday’s law of induction,

(2)

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,

(3)

where d_e is the electron inertial length, the corresponding terms in (3) taking into account finite-electron-mass 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 particle-in-cell (PIC) methods, but at the price of a large computational cost and memory requirements (when compared, for instance, to typical state-of-the-art hybrid-PIC 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 so-called 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 half-step 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 “advection-only advancement”¹ 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 second-order 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_eI (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:

(4)

where T_0, e is the constant, homogeneous electron temperature.

2.1.2. HVLF: Landau-fluid 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),

(5)

where bb is the dyadic tensor of the magnetic-field 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

(6a)

(6b)

where q_∥, e and q_⊥, e are the two scalars determining the parallel component of a gyrotropic heat-flux tensor (i.e., the flux of parallel and perpendicular thermal energy along the magnetic lines, respectively). In the hybrid-Vlasov code, these equations are used to advance in time p_∥, e and p_⊥, e employing a third-order Runge-Kutta method. Also, a sub-stepping 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 sub-steps, 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 full-kinetic 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 Chew-Goldberger-Low (CGL, Chew et al. 1956) double-adiabatic model they are set to zero (hence the name). Here, we adopt an LF type of closure for the electron heat-flux 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 three-pole and one-pole 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 heat-flux density coefficients, which in the case of the electrons are

(7a)

(7b)

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₀/m_pc is the proton cyclotron frequency (with B₀ the intensity of the magnetic field in the unperturbed state). In Eq. () the over-bar, , denotes the space average and ℋ is the operator represented, in Fourier space, by

(8)

where the hat, , denotes Fourier-transformed quantities. This operator effectively takes into account the distortion of magnetic-field 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 non-affordable computations in large simulations. Therefore, we adopted the semi-phenomenological 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 semi-phenomenological 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 ill-behaved 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 hard-coded “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 MHD-like models to perform numerical studies of magneto-rotational instability (Sharma et al. 2006), to perform numerical studies of the turbulent dynamo in the intracluster medium (Santos-Lima et al. 2014; St-Onge et al. 2020), and to derive a heating mechanism able to stably balance radiative cooling in the central regions of cold-core 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 electron-pressure components, LF closures (7a)–(7b) for the gyrotropic electron heat fluxes, and finite electron-inertia effect in the generalized Ohm’s law (3).

2.2. Full-kinetic model: iPIC

Compared to the LF closure implemented in the HVLF run, a full-kinetic model will provide also the non-gyrotropic effects of the dynamics of the electrons (i.e., those linked with FLR physics). Also, the quasi-neutrality 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 semi-implicit 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 semi-implicit method implemented on the code “iPic3D” (Markidis & Lapenta 2010; Lapenta 2012). The semi-implicitness 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 semi-implicit method has a reduced computational cost as compared with the implicit method and, clearly, with every explicit method. As a drawback, the semi-implicit 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 Vlasov-Maxwell 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 semi-implicit 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 quasi-neutrality 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 quasi-neutrality 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³ grid points. The iPIC (fully kinetic) run employs 14000 particle per cell (PPC), 7000 for each species. The hybrid simulations also include finite electron-inertia effects, and in all three models we adopted a reduced mass ratio of m_p/m_e = 100.

The initial setup consists of a “double-Harris (DH) sheet” configuration (Harris 1962) to accommodate for periodic boundary conditions:

(9)

where L₁ = 0.85d_p and L₂ = 1.7d_p are the two shear widths, y₁ = L₁L_y/[2(L₁ + L₂)] and y₂ = L_y − (L₂L_y)/[2(L₁ + L₂)] their positions, and B₀ = 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 non-vanishing 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 magnetic-field configuration in Eq. (9):

(10)

where the density profiles of the “background” and sheets populations, namely n_b, n₁, and n₂, respectively, are given by

(11)

and n₀ = 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₀L₁) and V_2, a = 2T_0, a/(B₀L₂) the velocities related to the current sheets. If one considers a single current sheet, L₂ → +∞, then Eqs. (9)–(11) form a consistent kinetic equilibrium (see Allanson et al. 2017). Due to the nonlinear nature of the Vlasov-Maxwell 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₁ ≃ n₀), the second one is negligible (n₂ ∼ 10⁻¹⁰n₀), and vice versa. Moreover, the choice of L₂ = 2L₁ makes the second current sheet quiescent over the entire simulation time (i.e., no reconnection nor other instabilities develop).

For our simulations, we considered²β_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 electron-FLR effects in the HVLF model. At the same time, the hybrid models employed in this work do include electron-inertia 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 magnetic-field 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 sub-steps), , 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 gyro-motion 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 well-resolved in time everywhere, with at least ≈30 time steps in the worst-case scenario (if τ_e = 2π/Ω_e is the gyro-period, 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 finite-k_⊥ρ_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 electron-Larmor-radius 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 pressure-anisotropy production and regulation.

In particular, we highlight the capability of the HVLF simulation of capturing part of the physics³ 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 O-point as follows (see Yeates & Hornig 2011):

(12)

The associated reconnection rate, R, is obtained as the time derivative of F:

(13)

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 time-averaging 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 “super-exponential” phase (up to ), and a final quasi-stationary 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, kt) (hence the name super-exponential or quasi-explosive, as defined in Ottaviani & Porcelli 1993, and Aydemir 1992), the third by a quasi-stationary 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 zoom-in 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 finite-PPC noise. Thus, we tested the convergence of the PIC simulation, in particular during the nonlinear stage, using different numbers of macro-particles, 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 late-time 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 “super-exponential” phase, while they begin to differ during the quasi-stationary regime. This feature is the consequence of the different response times of the electron dynamics associated with the three different models. During the quasi-steady 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 quasi-steady reconnection rate value around ∼0.04). The same trend for the quasi-stationary reconnection rates (i.e., higher with isotropic and isothermal electrons, lower in the full-kinetic 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 background-to-peak density ratio of n_b/n_H, 0 = 0.1 − 0.2⁴, in order to better compare with early simulations setup (usually GEM-like configurations, see Birn et al. 2001) and with some magnetospheric data. In Cassak et al. (2017) the so-called 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 background-to-peak density ratio is n_b/n₁(y₁) = 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 out-of-plane 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 quasi-stationary 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 ΔBz = Bz − Bz(t = 0) (first row), E∥ (second row), Jz (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 super-exponential phase. The plot is zoomed-in around the main X point. The is the rms value of E∥ in the shown region.

All the above-mentioned 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 double-adiabatic closure, with the addition of kinetic features⁵. This closure has been validated in Le et al. (2009, 2010b) and Egedal et al. (2013), within a two-fluid framework, and an analogous configuration for the out-of-plane 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 full-kinetic run. In this sense, the HVLF results are closer to those observed in the full-kinetic iPIC case.

Electron-anisotropy 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 steady-state 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 self-consistent 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 frozen-in 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

(14)

The first term on the right-hand side is the E × B-drift. The two other drifts are driven by the pressure-gradient 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 × B-drift 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 drift-current 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 quasi-stationary phase. The plot is zoomed-in around the main X point.

This resulting anisotropy-dominated 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⁶. In order to understand what balances this force, the steady-state electron momentum equation can be cast in the form

(15)

with F representing the other contributions (electric field, inertia, non-gyrotropy). Where the magnetic-field 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²| (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 B-curvature enhance electrons outflow jets, which generate the Hall magnetic field, B_H. This Hall component simultaneously provides the major contribution to the magnetic-field 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², 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 quasi-stationary phase⁷. 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 anisotropy-driven 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 proton-instability thresholds, as well as the electron PFHI curve, are computed for the maximum growth rate . The remaining electron-instability thresholds are computed for (which, in our simulations, means ).

Fig. 5. Cumulative histogram of the (β∥, a, T⊥, a/T∥, a) occurrences over the quasi-stationary 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.5dp. 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 fire-hose instabilities (green), and oblique fire-hose instabilities (violet). The black dash-dotted 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 full-kinetic case.

The main question remains why a non-negligible 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 time-cumulative 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 FHI-unstable 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 island-averaged points, which, in the HVM simulation, reach the lowest values among the three models. These features of the proton anisotropy-beta 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 quasi-stationary 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 island-averaged 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 MI-unstable.

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 magnetic-moment 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 anisotropy-related free-energy 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⁻²Ω_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 fastest-growing modes associated with these FHI branches cannot emerge due to the spatial constraints imposed by the islands size.

In Fig. 6, dash-dotted black lines are ∝1/β_∥, a curves, presented here to show how both electrons and protons populating the island seem to follow, on average, a double-adiabatic 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 quasi-stationary 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 full-kinetic 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 island-averaged point roughly follows a double-adiabatic 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 above-mentioned mirror-related 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 X-point 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 electron-FLR effects in the HVLF model). On top of that, we note that the portion of electron distribution populating the CI-unstable region in the HVLF simulation follows the 1/β_∥, e double-adiabatic 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 free-energy that would later lead to high-p_⊥,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 quasi-steady 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 anisotropy-limiting 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 quasi-stationary 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 magnetic-field 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 coarse-grained binning is employed for these plots.) The color-map refers to the bin-averaged 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 magnetic-field fluctuations are considered (i.e., ; Fig. 8, top row), one finds that the larger fluctuations are mainly located in the regions nearby marginal-stability 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 anisotropy-driven instabilities in their attempt to regulate electron-pressure 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, bin-averaged histogram of (top row), (middle row), and δB∥/δB⊥ (bottom row) in the (β∥, e, T⊥, e/T∥, e) plane. Histograms are cumulated over the quasi-stationary 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 so-called 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. Magnetic-fluctuations 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 non-propagating FHI, since at oblique propagation the fire-hose splits into a propagating, ω_r ≠ 0, and a non-propagating, ω_r = 0, branches (here, ω_r is the real-part 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 (non-compressional) mode, and it resonates with protons, not with electrons. The non-propagating 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 electron-mirror 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 magnetic-compressibility values in this region thus suggest that is the MI that tries to make up for the absence of the CI in regulating the electron-pressure 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 full-kinetic 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 electron-pressure anisotropy and anisotropy-driven (electron) instabilities, as well as the electron-kinetic response, such as Landau damping. Limitations that still hold within an HVLF approach with respect to a full-kinetic 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 so-called super-exponential phase. In the quasi-stationary stage, on the other hand, there is a quantitative difference in the quasi-steady reconnection rate between HVM and the other two models that is worth mentioning: The reconnection rate in the quasi-steady 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 quasi-steady reconnection rate is due to a background-to-peak 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 electron-pressure 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 anisotropy-versus-beta 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 parameter-space region that is beyond the OFHI marginal-stability 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 full-kinetic case, the anisotropic electron-fluid model with an LF closure efficiently regulates the distribution of electrons in an anisotropy-versus-beta 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 electron-LF model, which cannot capture the electron CI. The CI threshold is indeed effective at limiting the electron-pressure 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 electron-gyration 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 hybrid-kinetic framework (namely, including electron-pressure anisotropy and a fluid model of their linear kinetic response). When compared to the full-kinetic 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 full-kinetic 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 so-called transition range between MHD and sub-ion 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 electron-fluid model by including (large-scale) 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 full-kinetic response, although at a somewhat higher computational cost, and is currently considered to be a possible future direction.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 776262 (AIDA, www.aida-space.eu). F. F., S. S. C., and F. C. want to thank D. Telloni for useful discussions regarding in-situ 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 SuperMUC-NG 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 Max-Planck/Princeton Center for Plasma Physics (NSF grant PHY-1804048). 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 HPC-EUROPA3 (INFRAIA-2016-1-730897), 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

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 zoom-in 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 ΔBz = Bz − Bz(t = 0) (first row), E∥ (second row), Jz (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 super-exponential phase. The plot is zoomed-in 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 quasi-stationary phase. The plot is zoomed-in around the main X point.

In the text

Fig. 5. Cumulative histogram of the (β∥, a, T⊥, a/T∥, a) occurrences over the quasi-stationary 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.5dp. 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 fire-hose instabilities (green), and oblique fire-hose instabilities (violet). The black dash-dotted line is 1/β∥, a.

In the text

Fig. 6. Cumulative histogram of the (β∥, a, T⊥, a/T∥, a) occurrences over the quasi-stationary 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 island-averaged 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, bin-averaged histogram of (top row), (middle row), and δB∥/δB⊥ (bottom row) in the (β∥, e, T⊥, e/T∥, e) plane. Histograms are cumulated over the quasi-stationary 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