Influence of ion-to-electron temperature ratio on tearing instability and resulting subion-scale turbulence in a low-β e collisionless plasma

,


I. INTRODUCTION
Magnetic reconnection plays an important role in various space-plasma phenomena, from solar flares to geomagnetic storms.The investigation of magnetic reconnection has provided some crucial understanding of the mechanisms responsible for the release of energy particularly in the context of astrophysical plasmas, where the collisional mean free path is large enough for classical Coulomb collisions to be negligible.It is now acknowledged that reconnection in nature is often driven by collisionless effects.This necessitates models capable of including two-fluid effects, such as electron inertia.A significant step forward was taken in Refs.[2,38], where it was shown that, in the collisionless regime, two-fluid effects driving reconnection can provide a way to achieve fast reconnection.Subsequent bodies of work have confirmed the crucial role played by collisionless effects (see for example Refs.[5,7,25,28,32,37,53]), and have shown a good agreement with in situ spacecraft measurements in the Earth's magnetosphere [14].Notably, some of these studies involved fully kinetic simulations (as, for instance, in Ref. [22]).However, global fully kinetic simulations usually remain extremely expensive, and simplified models or hybrid approaches have emerged as alternatives with the potential to efficiently capture essential physical phenomena while significantly reducing computational costs.In the presence of a strong ambient magnetic field component, known as the 'guide field', gyrokinetic and gyrofluid models hold great potential for reconnection simulations, as shown in Refs.[16,37,46,51,54].
In this paper, we make use of a two-field gyrofluid model derived in Ref. [41] to simulate numerically reconnection events driven by electron inertia.This model isolates the dynamics of Alfvén waves (at the magnetohydrodynamics (MHD) scales) and kinetic Alfvén waves (at the sub-ion scales) in regimes where the couplings to the other kinds of waves are subdominant.It provides a good toolset for analyzing the plasma behavior in the strong guide field, low-β e regime, where β e is the ratio between electron kinetic pressure and guide field magnetic pressure.The model includes ion Larmor radius effects, and enables an arbitrary equilibrium ion-to-electron temperature ratio τ .For instance, this model is particularly valuable for investigating reconnection phenomena in regions such as the solar corona and its vicinity (β i = τ β e ≲ 0.1), the Earth's magnetosheath (where β i ≲ 1 can be observed), and the Earth's magnetosphere (β i ≪ 0.1) [52].Note that observations have revealed that, in most astrophysical plasmas, ion temperatures are usually larger than those of electrons, with for example τ ≳ 10 in the Earth's magnetosphere (see [11] for a measurement of magnetic reconnection at the magnetopause), τ ∼ 3 − 4 in the Earth's magnetosheath [21] and τ ∼ 2 in the solar wind [42].It thus appears relevant to investigate the effect of the temperature ratio within the framework of the two-field gyrofluid model, which is computationally less demanding than the kinetic descriptions.This model bridges the gap between reduced magnetohydrodynamics (RMHD), the inertial kinetic Alfvén waves (IKAW) model [14,39,40], and a reduced electron magnetohydrodynamics (REMHD) model that accounts for electron inertia, distinguishing it from the REMHD model derived in Ref. [48].In two spatial dimensions, the REMHD equations are formally identical to those of the electron magnetohydrodynamics (EMHD) model [6,34] which focuses on the incompressible regime, describing whistler waves.
The present work concentrates on the two-dimensional dynamics that develops in a plane perpendicular to the ambient field.Its aim is twofold.We first examine the linear growth rates of the tearing instability, investigating various equilibrium temperature ratios and confirming that the model converges toward the REMHD regime as τ increases.We point out that, compared to previous investigations on the role of the ion thermal radius and based on gyrofluid models [15,16,28], the present analysis does not require β i ≪ 1, where β i = τ β e .Relaxing this assumption makes it possible to access the above mentioned IKAW and REMHD regimes.We then study the turbulence regime resulting from reconnection in the cases of moderate and large values of the τ parameter.
Previous numerical simulations conducted in the cold-ion regime, using a reduced description that appears as a limit of our model, have provided evidence that collisionless magnetic reconnection can trigger fluid-like secondary instabilities [17][18][19][20]30].Always in the cold-ion limit, fluid-like secondary instabilities were observed in two and three-dimensional numerical simulations of a four-field model accounting for finite β e effects [31,33,50] It was observed that for low-β e reconnection, Kelvin-Helmholtz or Rayleigh-Taylor-like instabilities can develop, depending on the ratio of the ion-sound Larmor radius and electron skin depth ρ s /d e .
These secondary instabilities could potentially act as a source of turbulence.In the present study, we consider a broad range of values for the parameter τ and analyze the influence of this parameter on the nonlinear evolution of magnetic islands, as well as on the properties of the turbulence driven by the secondary instabilities.Nonlinear simulations were done for two distinct ion-to-electron temperature ratios, specifically τ = 100 and τ = 1.We consider these two cases representative of the finite-τ and large-τ regimes, respectively.In both regimes, we observed the existence of strong velocity shears that initiate Kelvin-Helmholtz instabilities.These instabilities lead to the propagation of turbulence through the separatrices and the formation of eddies.We will discuss the nature of this turbulence for both cases.
The paper is organized as follows.In Section II, we present the gyrofluid model and the different limiting regimes that it can cover.In Section III, we investigate the linear growth rates of the tearing instability in different parameter regimes.Section IV focuses on the turbulence dynamics and the vortex formation that develop at longer times.Section V is the Conclusion.

II. MODEL EQUATIONS
We make use of the gyrofluid model consisting of the two evolution equations complemented by the static relations 2 This system is the two-dimensional reduction of a model derived in Ref. [41].The latter is formulated in a slab geometry, adopting Cartesian coordinates {x, y, z} and assuming the presence of a strong magnetic guide field along the unit vector z.In the 2D version adopted here, we assume that the dynamical variables do not depend on the z coordinate.
Equations ( 1) and ( 2) correspond to the continuity equation for the electron gyrocenter density fluctuations N e and to Ohm's law, respectively.The relations (3) and ( 4), on the other hand, express the quasi-neutrality condition and the perpendicular component of Ampère's law, respectively.Here, A ∥ and B ∥ are the components, along the guide field, of the magnetic potential and of the magnetic fluctuations, whereas ϕ indicates the electrostatic potential.The expression for the total magnetic field is given by In Eqs. ( 1)-( 4) and ( 5), the variables are dimensionless and expressed according to the following normalization: where the hat denotes dimensional variables, L is a characteristic scale length, n 0 is the equilibrium density, B 0 the guide field amplitude, c the speed of light, the Alfvén speed, with m i indicating the ion mass, whereas di = c/ 4πn 0 e 2 /m i is the ion skin depth, with e corresponding to the elementary charge.This normalization differs from that adopted in Ref. [41].In the present paper, we opted for the normalization (6) because this will help in establishing contact with previous results present in the literature.Also, the normalization (6) might be more appropriate for astrophysical applications.
as a small expansion parameter, and assuming β e ∼ δ for δ → 0. We also recall that, although the low-β e limit tends to suppress electron FLR effects, our model retains one contribution which descends from an electron FLR term present in the parent gyrofluid model [9].This corresponds to the last term in Eq. ( 3), which becomes relevant in the large-τ limit, where it gets comparable to the retained contributions at scale d e .
We point out that, in the 2D version that we adopt here, the conservation laws differ considerably from those of the 3D version of the model.Indeed, as is typically the case with Hamiltonian reduced fluid models, in the 2D limit, the system acquires an infinity of Casimir invariants.For Eqs. ( 1)-( 4), these correspond to the two infinite families where C ± are arbitrary functions.All such functionals are conserved and constrain the dissipationless dynamics.The above infinite families include also quadratic functionals such as those used, for instance in Ref. [13], to investigate 2D drift-wave turbulence.Namely from a linear combination of quadratic Casimir invariants [41], one obtains in particular the generalized cross-helicity which is then also a conserved quantity for our 2D model.
In the general 3D case, on the other hand, instead of the infinite families C ± , one has the two linear Casimir invariants A generalized cross-helicity (with an expression analogous to Eq. ( 9)) is also conserved in 3D, but it is no longer a linear combination of Casimir invariants.
As already pointed out in Refs.[39,41], Eqs. ( 1)-( 4) generalize reduced models previously presented in the literature and which can be retrieved in the appropriate limits.In the following, we briefly review the different limits that will be relevant for the subsequent analysis.Specifically, under the conditions of finite k ⊥ ρ s and β e ≪ 1, the three regimes in τ are described in Subsections A, B, and C.These subsections correspond to scenarios where k ⊥ ρ i takes values of significantly less than 1, around 1, and greater than 1, respectively.
A. Cold-ion limit: τ ≪ 1 In this limit, the model ( 1)-( 4) reduces to where we have neglected contributions δ 2 times smaller than the leading order terms in each evolution equation, assuming τ = O(δ 2 ) and β e /2 ∼ d 2 e ∼ δ.In Eq. ( 13) we introduced the ), which is a modified ion sonic Larmor radius, accounting for parallel magnetic fluctuations effects.If β e is assumed even smaller, say β e = O(δ 2 ), then Eqs. ( 12)-( 13) identify with the two-field reduction, adopted for instance in Refs.[12,29], of the three-field model of Ref. [49].Note that, in this limit, in order to satisfy the relation (7), one also needs ρ s → 0 as δ → 0.

B. Finite τ and negligible β e
Here, contributions of order β e are assumed negligible in the model equations.Performing a Padé approximation of the operator Γ 0i (as done in Refs.[20,28]) one obtains, from Eq.
C. Hot-ion limit: τ ≫ 1 If we take τ ∼ 1/δ 2 , with β e /2 ∼ d 2 e ∼ δ and again neglect contributions O(δ 2 ) times smaller than the dominant terms in each evolution equation, we obtain which corresponds, up to the normalization, to the model for IKAW turbulence introduced in Ref. [14,40].
If we increase τ even further (say τ = O(1/δ 3 )), so that the contribution of 2/(τ β e ) becomes negligible with respect to 1, the evolution equations of the model reduce to Making use of Eq. ( 22) and introducing the whistler time t w = d i t, Eqs. ( 23)-( 24) can be rewritten as where b = d i B ∥ .Equations ( 25)-( 26) correspond to the equations of 2D REMHD, with electron inertia which, in two dimensions, are formally identical to those of 2D incompressible EMHD (see, e.g., Refs.[6,34]).We will therefore refer to the above limit as to the REMHD limit.However, later on, since the 2D equations are identical to those of EMHD, we will compare the numerical growth rates of the tearing modes to the theoretical formulas obtained within the framework of EMHD.We remark that, in the 3D case, the presence of the strong guide field ordering, would lead to a model for Inertial Whistler Turbulence [14].We also point out that the above procedure leading to REMHD upon neglecting corrections of order δ 2 makes it possible to retain corrections of order δ, such as those that become relevant at the scale d e .At the same time, it discards contributions at the scale of the electron Larmor radius ρ e = δρ s .Indeed, even if the only electron FLR term retained in the model leads to a contribution in Eq. ( 1) that is comparable to the other terms at scale d e when τ ≫ 1, it also leads in Eq. ( 2) to terms that become relevant at scales of order ρ e .Since other electron FLR contributions are missing in the model, the description at the scale ρ e is not accurate.
Therefore, it is necessary to ensure that the length scales involved in the solution remain significantly larger than the electron Larmor radius during the evolution.This issue will be further discussed in Sec.III B 2 in the context of tearing modes.

III. LINEAR TEARING MODES
In this Section, combining numerical and analytical studies, we analyze how the variation of the ion-to-electron temperature ratio τ affects the linear growth rate of tearing modes.
We consider both the so-called "constant-ψ" regime [26], typically valid for small values of ∆ ′ , and the regime with large ∆ ′ .We recall that the parameter ∆ ′ [26] is the variation, at the resonant surface, of the logarithmic derivative of the amplitude of outer solution in the linearized model equations and corresponds to the standard stability parameter of tearing modes.
In the following, we begin by describing the adopted equilibrium state, and subsequently treat the constant-ψ and large ∆ ′ cases.

A. Equilibrium state
We linearize the model equations ( 1)-( 2) about the equilibrium solution considered in Ref. [44] ϕ (0) = 0, A where ā is a constant determining the magnetic equilibrium amplitude.The corresponding poloidal equilibrium magnetic field is given by where y is the unit vector along the y coordinate.We set ā = 2.598, so that, for L x > 0.6586, one has max = 0 at equilibrium.The choice of the equilibrium profile is motivated by the rapid decay to zero of this profile and its derivatives for x → ±∞.This allows for the imposition of periodic boundary conditions, given a sufficiently large integration domain.This setting enables us to use a Fourier spectral code, most suitable for simulating gyrofluid models with nonlocal operators Γ 0 and Γ 1 , taking the form of Fourier multipliers, as demonstrated in [30].
This equilibrium corresponds to a current sheet centered at x = 0, with a dimensionless length of 2L y = 2 Ly /L, and a dimensionless width corresponding to unity.In other words, one uses the typical width of the current sheet as length unit.In our numerical simulations, the perturbation of the equilibrium potential ∥ (x, y, t), whose dependence on y is of the form A (1) ∥ ∝ cos(k y y), with k y = πm/L y (where m is an integer) and is initially excited by the mode m = 1.The stability condition is given by the tearing parameter, which for equilibrium ( 27) is (see [44]) The equilibrium ( 27) is unstable to tearing when ∆ ′ (k y ) > 0, corresponding to wave numbers In limiting cases, such as τ ≪ 1 and τ ≫ 1, where the gyroaverage operators become trivial, and the model equations become local in the space variables, other equilibrium profiles could be considered, such as the Harris sheet [5] or a sinusoidal profile [38].Such equilibria have a different expressions for ∆ ′ and therefore the corresponding linear growth rates of the tearing instability display a different dependence on k y .
The linear and maximum growth rate of the tearing instability are measured at the X-point using the following quantity: which is constant in time during the linear phase.

B. "Constant-ψ" regime
In the "constant-ψ" regime, it is assumed that the amplitude of the perturbation A (1) ∥ is approximately constant in the inner region centered around the resonant surface.This assumption is valid when ∆ ′ ϵ ≪ 1, where ϵ denotes the width of the inner region, which is estimated below.
We first present analytical results concerning the cold-ion and EMHD regimes, that will then be compared with the numerical values obtained by the simulations in the appropriate limits.

Dispersion relation for the cold-ion regime
We consider Eqs. ( 12)-( 13) and their linearization about the equilibrium (27).We assume perturbations of the form (x, y, t) = φ(x) exp(ik y y + γt) + c.c., (31) where c.c. indicates complex conjugate.In Eq. ( 31), the constant γ provides the growth rate of the tearing perturbation.We assume that the amplitudes Ã and φ are even and odd functions of x, respectively, which is the standard parity for tearing modes.
When β e is neglected compared to unity, it is known that the dispersion relation for tearing modes for the considered equilibrium is given by [4,43] with ∆ ′ given by Eq. ( 29).
We remark that small corrections in β e due to parallel magnetic perturbations only affect Eq. ( 13), in particular by turning ρ 2 s into ρ ′ s 2 .Therefore, such effects can be accounted for, in the dispersion relation, in a straightforward way, by replacing ρ s with ρ ′ s in Eq. (32).
Because ρ ′ s ≤ ρ s , we can infer that parallel magnetic perturbations tend to decrease the growth rate of tearing modes.We recall that a similar result was obtained in Ref. [27].In this Reference, the combined effects of parallel magnetic perturbations and electron FLR terms were considered, but the latter ones were not treated self-consistently, although the corresponding relation showed good agreement with numerical results.The present model neglects electron FLR terms in the cold-ion limit but retains small corrections due to parallel magnetic perturbations.Therefore, it makes it possible to isolate and identify, in a more rigorous way, the stabilizing effect of a finite but small B ∥ .

Dispersion relation for the EMHD regime
For the EMHD model in the constant-ψ regime, an analytical dispersion relation for tearing modes is given in Ref. [10].A thorough investigation of EMHD tearing modes, including a numerical validation of the dispersion relation of Ref. [10], was recently presented in Ref. [3].However, applying the result of Ref. [10] to our gyrofluid model in the large τ limit is more delicate.As anticipated in Sec.II C, the cold-ion and EMHD systems, corresponding to Eqs. ( 12)-( 13) and (up to the normalization) ( 23)-( 24), respectively, were obtained from the original gyrofluid model, assuming that length scales remain of the order of the characteristic scale length L. The latter, recalling that x = x/L and considering Eq. ( 28), can be considered as a characteristic scale of variation of the equilibrium magnetic field.
However, in the tearing mode analysis, an inner region, involving small scales, is present around the resonant surface, where terms with high order derivatives, that are negligible on scales of order L, can become relevant.In the cold-ion case, one can safely insert the relations ( 14)-( 15) into Eqs.( 1)-( 2), neglect corrections which are δ 2 smaller and thus obtain Eqs. ( 12)-( 13), for which the dispersion relation ( 32) is valid.Indeed, it thus turns out that terms which are δ 2 smaller than the leading terms, are negligible both in the inner region and in the remaining outer region, without imposing further conditions between the growth rate and the parameters β e , d e and ρ s .However, for the hot-ion case, in general, this is not the case.Therefore, we first examine the conditions under which the EMHD tearing mode relation of Ref. [10] can be applied to the hot-ion limit of our gyrofluid model.
Because we are interested in the EMHD limit, rather than in the IKAW limit, we assume that τ is large enough, so that 2/(τ β e ) is neglected compared to 1 in Eq. ( 21).Then we insert the resulting hot-ion relations into Eqs.( 1)-( 2) but without neglecting any terms, so as to account for the possibility that terms negligible on scales of order L might become non-negligible in the inner region.Upon linearizing the resulting system about the equilibrium (27), we obtain where the prime denotes derivative with respect to the argument of the function, which is x in this case.We consider the system ( 35)-( 36) on an infinite domain (which is a reasonable approximation if L x is sufficiently large) with Ã → 0 and φ → 0 for x → ±∞.
The EMHD case is retrieved from Eqs. ( 35)-( 36) if the last three terms on the left-hand side of Eq. ( 36) are neglected and the change of variables γ = d i γ w , φ = d i b is performed.
In the latter expressions we indicated with γ w the growth rate in terms of whistler units of time and with b the amplitude of the perturbed parallel magnetic field.
Assuming γ ≪ 1 and d e ≪ 1, in the outer region, far from the resonant surface at x = 0, all terms proportional to d 2 e and γ are negligible.These include also the three above mentioned "non-EMHD" terms, which confirms that, on scales of the order of the equilibrium magnetic field, the tearing analysis of the system (35)-( 36) coincides with that of EMHD.However, around x = 0, where A where is a stretching factor and Given that ϵ ≪ 1, around x = 0 the system ( 35)-( 36) can be approximated by where we made use of the constant-ψ approximation and set Ãin ≈ 1.Two out of the three non-EMHD terms turned out to be negligible also in the inner region.However, a term remains, which is given by the last term on the left-hand side of Eq. ( 41).This term originates from the electron FLR effect retained in Eq. ( 3) and is negligible, compared to the third term on the left-hand side of Eq. ( 41), if where ρ e is the thermal electron Larmor radius.Therefore, the EMHD case is retrieved if the width of the inner region is much larger than the electron Larmor radius.
Note that, in the cold-ion case, retaining the electron FLR correction, from Eq. ( 3), leads to N e = (1 + δ 2 )∇ 2 ⊥ ϕ.Therefore, the electron FLR term, in this case, corresponds to a small correction of order δ 2 to the ion polarization term, and can thus be safely neglected, which leads to the relation (14).However, the requirement that the width of the inner region should be much larger than the thermal electron Larmor radius, should hold also in the cold-ion limit, as pointed out in Refs.[24,27].
The condition (42) can be made more explicit once one finds the expression for ϵ in terms of the parameters of the system.This corresponds to the distinguished limit of the system.
The second term on the left-hand side of Eq. ( 43) would be of the same order of the first term, if ϵ were equal to d e .However, this would make the third term diverge, because its coefficient would become d 2 i ā2 X 2 /g 2 , which goes to infinity as g → 0. Consequently, the second term has to be subdominant.We remark that, if we had kept the contribution proportional to 2/(τ β e ) of the hot-ion limit, the second term on the left-hand side of Eq. ( 43) would have been multiplied by a factor (1 + 2/(τ β e )).If this factor remains of order unity, that term would have still been subdominant.Therefore, the presence of the coefficient characteristic of the IKAW regime does not affect the tearing dispersion relation, which is consistent with the analysis of Ref. [8].
The distinguished limit is then provided by balancing the first and the third term, which leads to The condition (42) can then be reformulated, in terms of the original parameters of the system, as If the condition ( 45) is fulfilled, the tearing analysis of Ref. [10] applies to the hot-ion limit of our gyrofluid model, which yields where Γ(x) is the Gamma function.We also remark that, with respect to the original Ref. [10], the expression (46) contains a factor d i , due to the fact that γ is normalized in terms of Alfvén units of time.
Assuming the validity of the relation (46) and inserting it into Eq.( 45), one can also re-express the latter condition as In this formulation, the condition for excluding the electron FLR contribution only depends on the plasma parameters β e and d e and on the parameters ā and ∆ ′ associated with the equilibrium and the wave number of its perturbation.

Numerical results and comparison with analytical dispersion relations
Here we analyze, by means of numerical simulations, the tearing growth rates for different values of the parameters and compare the results, in the appropriate limits, with the analytical predictions of Secs.III B 1 and III B 2. Equations ( 1) -( 4) are solved with a Fourier pseudo-spectral solver on a periodic domain, with advancement in time performed by a third-order Runge-Kutta numerical scheme.To investigate the linear growth rate, we consider the case of small ∆ ′ , where the "constant-ψ" approximation is applicable.In this context, the initial wave number is set to k y = 1.923, leading to ∆ ′ = 1.699.
To study the influence of τ , we conducted a series of numerical simulations covering a range of equilibrium temperature ratios spanning from τ = 0.1 to τ = 1000.In terms of the ion thermal Larmor radius ρ i = √ τ ρ s , the range corresponds to going from ρ i = 0.09 to ρ i = 9.48.For these simulations we fixed the sound Larmor radius to ρ s = 0.3 and the mass ratio to m e /m i = 2 × 10 −5 .The use of such a small mass ratio, which is made possible by the computational efficiency of the gyrofluid modeling, ensures a small electron skin depth and a small β e , meeting the validity criteria for the linear theory.
Figure 1 displays the linear growth rates for three different values of β e and d e given by simulations together with the asymptotic values predicted by the theory.The range of parameters has been chosen in such a way that, for τ ≫ 1, condition ( 45) is satisfied, so that the EMHD regime should be recovered.For τ = 0.1, the numerical growth rates closely align with those predicted by the cold-ion relation (32), shown in red on the figure (the adopted values of β e in this case are so low that the use of the effective parameter ρ ′ s , mentioned in Sec.III B 1, is practically of no use).As we increase τ , the growth rates also increase and the model converges towards the EMHD model described in Sec.II C. For τ = 1000, the linear growth rates match the analytical estimates based on Eq. ( 46) and represented in blue on the figure.
Although the values of τ of practical physical interest, for instance for space plasmas, concern only a small portion around τ ≈ 1 of the considered range, our analysis shows how,  46) and corresponding to the EMHD regime.
by increasing τ , the linearized system can perform a transition from the cold-ion regime, to a regime where the thermal ion Larmor radius is much larger than the characteristic scale of the current sheet, so that the observed dynamics is essentially due to electrons.In this regime the predictions of EMHD apply.In particular, we note the transition from a growth rate linear in d e ∆ ′ , for cold-ions, to a growth rate proportional to d i d 2 e ∆ ′ 2 for very hot ions.
In Tables I and II, we present a precise quantitative comparison indicating the errors between numerical and theoretical values.Table I shows the results obtained for τ = 1000 compared to the EMHD theoretical formula (46), while Table II shows the numerical results for τ = 0.1 compared to the theoretical predictions based on Eq. ( 32).The values of the growth rate displayed in the Tables correspond mostly to those shown also in Fig. 1.The applicability conditions of relations ( 32) and (46) appear to considerably limit the admissible range of the plasma parameters, if one requires a relative discrepancy of a few percent between numerical and analytical predictions.However, we also observe a satisfactory agreement with the EMHD relation (46) in the case d e = 0.21 and β e = 2 × 10 −5 , well different from those shown in Fig. 1.
We also verified how the agreement between the numerical growth rate and the EMHD growth rate obtained from Eq. ( 46) improves as the condition ( 45) is better and better fulfilled in an asymptotic sense.Figure 2 provides the evolution of γ as a function of 1/β e .
In this case, we maintain τ = 1000, and d e = 0.21, and we vary β e and ρ s , while keeping d e and d i constant.If one assumes the expression (46) for the growth rate, the condition (45 Given that we keep d e and ∆ ′ fixed, this condition is better and better fulfilled as β e de- Num.value Theo.EMHD value FIG.2: Linear growth rate as a function of 1/β e for τ = 1000.As β e decreases, the numerical growth rate approaches the EMHD analytical prediction (46) (solid line), also reproduced by direct simulation of the EMHD equations (x-cross).This reflects the fact that one is approaching the asymptotic regime of the condition (48).

C. Large ∆ ′ regime
We now focus on simulations initiated with a wavenumber of k y = 0.55, yielding ∆ ′ = 48.45.The ion-sound Larmor radius is fixed to ρ s = 1.On Fig. 3 we report the linear growth rate as a function of τ for both β e = 0.005 and β e = 0.02.
The nonlinear evolution of two of these runs, both with β e = 0.02, is discussed in detail in the next section.
The evolution of the growth rate follows a trend similar to that observed in the case of small ∆ ′ , with γ increasing as τ increases.Nevertheless, some differences can be identified.
For τ = 0.1, the scaling with respect to d e is rather close to the scaling γ ∝ d 1/3 e predicted in Ref. [43].Indeed, the ratio between the growth rates for d e = 0.44 and d e = 0.22 is 0.77/0.63≈ 1.222, while (0.44/0.22) 1/3 ≈ 1.259.A direct comparison with analytical predictions can be made using the formula [43] which is valid also for large ∆ ′ .
The results of the comparison are summarized in Table III.One can notice that the  relative errors are greater than those in Table II.In fact, whereas in Sec.III B the values of the parameters were chosen in order to approach the asymptotic regime of validity of the theory, for the large ∆ ′ case, the emphasis is mainly in complementing the nonlinear analysis of Sec.IV, for which the choice of the parameters was not dictated by adherence to analytical linear theory.Therefore, a greater discrepancy between numerical and analytical results could be expected.In particular, for the case d e = 0.44, we see that, for instance, the condition γ ≪ 1 is not well fulfilled.
With regard to the hot-ion case, i.e. τ = 100, we could argue that it falls into the EMHD regime, similarly to what we did for the "constant-ψ" case.The most updated tearing linear theory against which we could test this, is provided, to the best of our knowledge, in Ref.
[1] and predicts, using Alfvén units of time, γ e .However, the ratio between the two growth rates of Fig. 3, for τ = 100, is 6.6/3.61≈ 1.828, whereas, according to the theory of Ref. [1], their ratio should be (20/10) × 2 2/3 ≈ 3.175.Thus, we are evidently out of the regime of validity of the theory of Ref. [1] and what we observe is a weaker increase of the growth rate.We remark that in Ref. [3] a disagreement between such theory and numerical results of direct EMHD simulations was found, suggesting that this theory is quite sensitive to the adopted regime of parameters.On the other hand, numerical results of EMHD simulations in Ref. [1] show, in some regimes, scalings other than those predicted by the theory, and which could be more compatible with what we observe for τ = 100.
In summary, in the regime of small ∆ ′ , when increasing τ , we observe a compelling convergence of the numerical growth rate to the EMHD growth rate obtained from Eq.( 46), provided condition (45) be satisfied.Figure 2 also demonstrates a remarkable convergence, as β e decreases, of the numerical growth rate towards the asymptotic prediction.Differently, in the regime of large ∆ ′ with hot ions, a noticeable deviation from the predictions of the theory of Ref. [1] is observed, indicating a strong sensitivity to parameter regimes, a phenomenon already mentioned in other studies [3].

IV. TURBULENCE GENERATED BY RECONNECTION
In this Section we discuss regimes in which Kelvin-Helmholtz instabilities, following an initial tearing instability, lead to a turbulent current layer.The simulations considered for this study were conducted using a grid size of 2080 2 collocation points in a 2D domain defined as −2.2π ≤ x ≤ 2.2π and −1.8π ≤ y ≤ 1.8π.For these runs, the characteristic length L is taken equal to ρs .We applied filters [35] designed to smooth out scales for which k⊥ L = k ⊥ > 300.This range is well beyond the region in which k ⊥ ρ e ∼ 1, suggesting that the filter has a minimal impact on the spectral domain under investigation.Our primary interest concerns the dynamics at the electron inertial scale d e .
The two simulations presented in this Section were performed with the following param- Similar secondary fluid-like instabilities taking place at the separatrix of a magnetic island were also observed in the context of tearing Particle-In-Cell simulations [23,45].In the cold-ion case, a secondary instability taking place at the separatrices was reported in Ref. [47].Fig. 5 shows the out-of-plane electron gyrocenter velocity, revealing the presence of large and small-scale current structures.For comparison, we superimpose in the right panel, a blue square with side of length πρ e equal to the half wavelength corresponding to the wavenumber k ⊥ = 1/ρ e .The size of the vortices is considerably greater than the size of the square.Therefore it is reasonable to conclude that the formation of such vortices does not FIG.6: Same as in Fig. 5 for the REMHD simulation.
depend on the dynamics occurring at the scale of the electron Larmor radius, where our model is not accurate.It is of interest to compare this run with an integration of the 2D REMHD model ( 25)-( 26).The plots showing the time evolution of the current are presented in Fig. 6.We note a qualitative similarity with the simulation with τ = 100, in particular the development of a secondary Kelvin-Helmholtz instability and the subsequent formation of magnetic vortices.It can be observed that the magnetic vortices (dashed black lines, corresponding to the isolines of A ∥ ) persist and have a similar shape in the two simulations.
Figure 7 presents cuts through a plasmoid in the parallel electron velocity for the two runs.
These structures have been fitted, for comparison, with a Gaussian whose full width at half maximum is FWHM = 0.3 in units of L = ρs for both τ = 100 and REMHD simulations, indicating that the size of the vortices is typically 1.3 de .We see the onset of a "double structure" especially conspicuous in the REMHD, with the internal structure becoming sharper when the dissipation is reduced and possibly singular in the zero-dissipation limit (not shown).This issue deserves further investigations.close to −4.The slope of −4 is not in significant disagreement with the −11/3 ≈ −3.67 spectral exponent predicted for IKAW (inertial kinetic Alfvén wave) turbulence [14,36], but the presence of structures formed by the Kelvin-Helmholtz instability can certainly impact the perpendicular magnetic spectrum slope.We note that, the spectrum in the region for k ⊥ < 10 may be subtly influenced by the presence of the equilibrium magnetic field, while the knee appearing around k ⊥ = 300 corresponds to the impact of filters.Secondly, A detailed study of this issue would preferably be performed in the framework of a homogeneous isotropic turbulence.In the case with τ = 1 (corresponding to an equilibrium current sheet with a width of the order of the ion Larmor radius), Fig. 9 shows that the velocity shear in U ⊥ is less important compared to the case τ = 100 (which corresponds to a current sheet of a width 10 times smaller than the ion Larmor radius) leading the Kelvin-Helmholtz instability to develop at a smaller scale.When considering the electron gyrocenter velocity (Fig. 5), we observe we conducted a simulation in which we deliberately removed the electron FLR term from the code (note that k ⊥ ρ i remains unchanged in this comparison, unlike in the previous case of hot ions).The resulting out-of-plane electron gyrocenter velocity is presented in Fig. 11, highlighting a notable reduction of magnetic vortices and an absence of current structures of sizes similar to d e or smaller.On the other hand, for this run, we had to increase the range of the spectrum affected by the filter, resulting in smoother smallest scales.We can see that, for this run, the region k ⊥ > 100 is affected by the filter.However, given the high resolution of our simulation, there remains a sizeable range below the scale ρ e .This observation could suggest the existence of a physically significant dynamics at the scale of the electron Larmor radius, which lies beyond the scope of the present model.Retaining all electron FLR terms could be of interest for a future work, but, at finite ion temperature,

V. CONCLUSIONS
In this study, we have investigated the linear and nonlinear evolution of the tearing mode instability in a collisionless plasma for different ion-to-electron temperature ratios.Through numerical simulations, we explored the behavior of the instability in different regimes, focusing on cases with moderate and large ion Larmor radii compared to the characteristic width of the equilibrium current sheet.In particular we studied the transition of the instability behavior from a cold-ion to an EMHD regime as τ increases.In the small ∆ ′ regime, we numerically validated analytical dispersion relations and showed how the growth rate evolves from a linear scaling in d e ∆ ′ for cold ions [44] to a scaling in d 2 e ∆ ′ 2 for hot ions [10].
As expected, the agreement with the analytical EMHD growth rate requires the width of the inner region be greater than the electron thermal Larmor radius.
In the regimes τ = 100 and τ = 1, we examined the development of turbulence triggered by the Kelvin-Helmholtz instabilities in the nonlinear phase of the tearing mode, concentrating on the influence of the ion-to-electron temperature ratio on the development of turbulence.
For τ = 100, the presence of large magnetic vortices, formed as a consequence of the Kelvin-Helmholtz instability, is noted.The size of such vortices is much larger than the scale at which electron FLR effects (which are not consistently taken into account in our model) become relevant.We can then assume that their formation is genuine and not influenced by deficiencies of the model.We observed analogous vortices also in EMHD simulations, where electron FLR effects are absent.At a closer inspection, EMHD vortices exhibit a "double structure", with an inner core whose width appears to be influenced by dissipative effects.The turbulent spectrum exhibits power law of −4, which requires further investigation, preferably in the framework of homogeneous and isotropic turbulence.
For τ = 1, and β e = 0.028 we also observe a secondary Kelvin-Helmholtz-like instability with subsequent formation of vortices.However, in this case, the electron FLR term present in our model appears to play a more important role in their formation.Concerning the spectrum, a power law close to −13/3 is observed.In this regime, the presence of small-scale structures, their reduction upon the removal of the electron FLR term, and the unexplored dynamics at scales below ρ e emphasize that further explorations, with the inclusion of all the electron FLR effects and the parallel ion dynamics, could unlock valuable insights.
Future developments of this study also include simulations in three space dimensions, where a richer and more complex dynamics is expected.Also in the perspective of including the dynamics of higher order moments and kinetic effects, the 2D assumption can be rather restrictive.For instance, as shown in Ref. [51], kinetic effects on the tearing growth rate, described by means of a Landau closure, become appreciable only in 3D.Furthermore, investigating reconnection arising during the evolution of a freely-decaying turbulence would allow us to explore the self-organization of the plasma in the absence of external driving.
These extensions will contribute to a deeper understanding of the interplay between reconnection and turbulence in magnetized plasmas.

s k 2 ⊥ , with k 2 ⊥
Parameters of the system are the electron skin depth d e = (1/L)c/ 4πn 0 e 2 /m e , with m e indicating the electron mass, the ratio β e = 8πn 0 T e0 /B 2 0 between electron kinetic and guide-field magnetic pressures, the ion sonic Larmor radius ρ s = β e /2 m i /m e d e and the ion-to-electron equilibrium temperature ratio τ = T i0 /T e0 .The model equations also involve the canonical bracket [f, g] = ∂ x f ∂ y g − ∂ y f ∂ x g and of the ion FLR operators Γ ni (b i ), for n = 0, 1, which correspond, in Fourier space, to multiplication by I n (b i ) exp(−b i ), where I n is the modified Bessel function of first type of order n and b i = τ ρ 2 indicating the squared modulus of the wave number in the plane perpendicular to the guide field.We also indicated with ∇ 2 ⊥ the perpendicular Laplacian operator defined by ∇ 2 ⊥ f = ∂ xx f + ∂ yy f .The independent variables of the model vary on a domain D = {(x, y) : −L x ≤ x ≤ L x , −L y ≤ y ≤ L y }, where L x and L y are positive constants.Periodic conditions are imposed at the boundaries of D.This model can describe the dynamics of collisionless plasmas at low-β e , accounting for ion FLR effects as well as electron inertia, which can break the frozen-in condition and allow for magnetic reconnection.The electron fluid is assumed to be isothermal, whereas the fluctuations of the ion gyrocenter moments are neglected in Eqs.(1)-(4) −Lx≤x≤Lx|B eq⊥ (x)| = 1.Note also that, from Eqs. (3)-(4), the choice ϕ(0) negligible in the outer region might become important.To see at which scale such terms become non-negligible one introduces rescaled variables

FIG. 1 :
FIG. 1: Linear growth rate as a function of τ for different values of β e and d e .The normalized ion Larmor radius becomes larger than unity starting from τ = 10.Red symbols indicate the theoretical predictions based on the relation (32) valid in the cold-ion case.Blue symbols represent the theoretical predictions based on Eq. (46) and corresponding to the EMHD regime.
creases.The plot includes a theoretical reference line representing the growth rate value derived from the EMHD theory.Interestingly, as β e decreases, corresponding to smaller ρ s values, our model exhibits a closer agreement with the EMHD results.The deviations from the theoretical EMHD value as β e increases clearly demonstrate the significant impact of β e .

22 FIG. 3 :
FIG. 3: Linear growth rate as a function of τ for different values of β e and d e , when ∆ ′ = 48.45.

FIG. 5 :
FIG. 5: Color-scale plot of U e with isolines of A ∥ showing the evolution of the turbulence outside the separatrix for τ = 100.The displayed region corresponds to approximately a quarter of the total domain, focusing on the separatrix of the island.The left panel also includes a zoomed-in view within the square area delimited by the grey line, where fluid vortices start developing.In the rightmost panel a square of size πρ e is superimposed.This shows that the size of the vortices is considerably greater than the scale at which electron FLR effects become relevant.

Time-averaged spectraFIG. 7 :
FIG. 7: Cut of U e along x for the run τ = 100 (left panel) and for the REMHD run (right panel).Both structures are fitted by a Gaussian (orange dashed lines).

FIG. 8 :
FIG. 8: Left: Time averaged spectra of the field B ⊥ for τ = 100 and for the REMHD run.

FIG. 9 :
FIG. 9: Left panel: Color-scale plot of U y for τ = 1.Right panel: Profiles of U ⊥ and B ⊥ along the black arrow visible in the right panel.

TABLE I :
Comparison, for τ = 1000, of the growth rate γ obtained numerically and γ T , based on Eq. (46), for different values of the parameters.Relative errors are also indicated.

TABLE II :
Comparison, for τ = 0.1, of the growth rate γ obtained numerically and γ T , based on Eq. (32), for different values of the parameters.Relative errors are also indicated.

TABLE III :
Comparison, for τ = 0.1, of the growth rate γ obtained numerically and γ T based on Eq. (49), for different values of the parameters.Relative errors are also indicated.