E Tassi

M Romé

C Chandre

A Hamiltonian fluid-kinetic model for a two-species non-neutral plasma

A model for describing the dynamics of a pure electron plasma in the presence of a population of massive charged particles is presented. The model couples the fluid dynamics of the pure electron plasma with the dynamics of the massive particle population, the latter being treated kinetically. The model is shown to possess a noncanonical Hamiltonian structure and to preserve invariants analogous to those of the two-dimensional (2D) Euler equation for an incompressible inviscid fluid, and of the Vlasov equation. The Hamiltonian structure of the model is used to derive a set of stability conditions for rotating coherent structures of the two-species system, in the case of negatively charged massive particles. According to these conditions, stability is attained if both the equilibrium distribution function of the kinetic species and the equilibrium density of the electron fluid are monotonically decreasing functions of the corresponding single-particle energies in the rotating frame. For radially confined equilibria near the axis, the stability condition corresponds to the existence of a finite interval of rotation frequencies for the reference frame, with the upper bound determined by the presence of the kinetic population.

The dynamics of single-species non-neutral plasmas has been studied extensively in the literature [START_REF] Davidson | An Introduction to the Physics of Nonneutral Plasmas[END_REF][START_REF] O'neil | [END_REF]. We consider here the behavior of a multi-species non-neutral plasma, and more specifically an electron plasma, subject to a uniform and constant magnetic field, and in the presence of a population of massive negatively charged particles.

An example of a multispecies negatively charged nonneutral plasma is the electron-antiproton plasma used as an intermediate step in the production of trapped antihydrogen [3]. The experimental investigation of the dynamics of a non-neutral complex plasma, and in particular a plasma of electrons contaminated with a small fraction of sub-micrometric charged dust particles is the subject of the "DuEl" project [START_REF] Romé | Non-Neutral Plasma Physics VIII[END_REF][START_REF] Maggiore | [END_REF].

Under conditions of strong magnetization of the electron component, which are typical of the experiments performed on non-neutral plasmas confined in Penning-Malmberg traps, a cold fluid description of the electron dynamics is adopted, in which electron inertia as well as temperature effects are neglected [START_REF] Davidson | An Introduction to the Physics of Nonneutral Plasmas[END_REF]. On the other hand, the dynamics of massive charged particles occurs on slower time scales, since their charge to mass ratio is orders of magnitude smaller than the corresponding ratio of the electrons. Moreover, at relatively low magnetic field strengths regimes can be found where ions are only weakly magnetized.

Under these circumstances, the simplest model to describe these multi-species systems consists of a hybrid system which couples, through Poisson's equation, a kinetic equation for the massive species with a fluid equation for the electron guiding centers. A further simplification that can be adopted for these systems consists of considering a 2D dynamics. For pure electron plasmas in Penning-Malmberg traps the high-frequency longitudinal electron motion can be averaged out for a wide range of operational regimes, and the electron plasma dynamics becomes isomorphic to that of a 2D incompressible inviscid fluid, where the fluid vorticity corresponds to the electron density and the fluid stream function to the electrostatic potential [6,7]. A 2D multi-fluid model turns out to be appropriate also for a cold non-relativistic electronion system at high magnetic field strengths [START_REF] Davidson | An Introduction to the Physics of Nonneutral Plasmas[END_REF], while a kinetic treatment for the dynamics of the ion population is required at moderate values of the confining field. The hybrid model we present here can obviously be generalized to account for multiple negative ion species treated kinetically. Furthermore, the model can also be applied to the case of a positron plasma in the presence of kinetic populations of positive ions (see, e.g., Ref. [8]).

On the basis of the above assumptions we consider a plasma contained in a cylindrical domain of radius a and length L z and we use cylindrical coordinates (r, θ, z), with z along the constant and uniform magnetic field, B = B 0 ẑ. Our model equations are then given by ∂n e ∂t + [φ, n e ] = 0, (1)

∂f ∂t + v • ∇f + q I M -∇φ + v c × B + g • ∂f ∂v = 0, (2) ∆φ = 4π en e -q I d 2 vf . (3)
In the above system n e indicates the guiding center electron density, f (r, θ, v r , v θ ) is the distribution function of the massive particle population (integrated over v z ), φ is the electrostatic potential, -e the electron charge, q I < 0 the charge of the massive particle, M its mass, g = -g(sin θr + cos θ θ) the acceleration due to gravity (the latter may become non-negligible for heavy ions), c the speed of light and d 2 v = dv r dv θ the surface element in velocity space. Finally, we defined the bracket [f, h] between two functions f and h of the spatial variables (r, θ)

as [f, h] = (c/B 0 r)(∂ r f ∂ θ h -∂ θ f ∂ r h).
Equation ( 1) represents a continuity equation for the electron guiding center density, which is advected by the incompressible E × B flow. Equation ( 2) is a Vlasov equation that governs the evolution of the massive particle population, whereas Eq. ( 3) is the Poisson's equation determining the electrostatic potential.

Equations ( 1)-( 3) form a Hamiltonian system assuming that boundary terms vanish when integrating by parts. The Hamiltonian functional H and the Poisson bracket {, } are given by

H(n e , f ) = dAd 2 v M v 2 r + v 2 θ 2 + M gr sin θ f - 1 8π dAφ∆φ, (4)
and

{F, G} = dA n e e [F ne , G ne ] + dAd 2 vf [F f , G f ] v , (5)
respectively. In Eq. ( 5), the subscripts of F and G indicate functional derivatives, whereas dA = rdrdθ is the surface element. The bracket [f, h] v in Eq. ( 5) is defined by

[f, h] v = 1 M ∂f ∂r ∂h ∂v r - ∂f ∂v r ∂h ∂r + 1 r ∂f ∂θ ∂h ∂v θ - 1 r ∂f ∂v θ ∂h ∂θ + v θ r + q I B 0 M c ∂f ∂v r ∂h ∂v θ - ∂f ∂v θ ∂h ∂v r . (6)
Note that in Eq. ( 4) the Hamiltonian H is defined as an energy per unit length. It can be easily verified that the dynamical equations ( 1) and ( 2) are obtained from the expression ∂χ/∂t = {χ, H} which governs the evolution of a generic Hamiltonian system, by replacing χ with n e and f , the Hamiltonian functional with Eq. ( 4), the Poisson bracket with Eq. ( 5), and making use of the Poisson's equation (3). Note that the bracket in Eq. ( 5) is the direct sum of two independent Poisson brackets, namely those for the incompressible 2D Euler equation [START_REF] Morrison | Mathematical Methods in Hydrodynamics and Integrability in Related Dynamical Systems[END_REF] and for the Vlasov equation [10,11]. Consequently it qualifies as a Poisson bracket, and in particular it satisfies the Jacobi identity. This bracket possesses two independent families of Casimir invariants, that is functionals F , such that {F, C} = 0 for all functionals C(n e , f ). These families are given by

C 1 = dAF(n e ), C 2 = dAd 2 vG(f ), (7)
with F and G arbitrary functions. The invariants include, as particular cases, the familiar invariants of the 2D Euler incompressible equations, such as the enstrophy, and those of the Vlasov equation, such as the entropy and the total number of particles. We remark that further examples of Hamiltonian hybrid fluid-kinetic models, accounting for coupling via the current density, or via the pressure, have been discussed in Refs. [START_REF] Morrison | Nonlinear Physical Systems Spectral Analysis, Stability and Bifurcation[END_REF][START_REF] Tronci | [END_REF][14].

The Hamiltonian structure of the model is exploited to obtain stability conditions for equilibria of the hybrid system (1)-(3) using the Energy-Casimir method (see, e.g., Refs. [15][16][17][18]). We restrict the analysis to the case of two species rotating azimuthally with the same uniform and constant frequency ω. Equilibria of interest are then those which are stationary in the rotating frame. In order to apply the Energy-Casimir method, the Hamiltonian structure of the model equations in the rotating frame is required. A convenient way to accomplish this is to look for a new set of field variables, in terms of which the Poisson bracket keeps the same form, whereas additional contributions accounting for the change in the reference frame appear in the Hamiltonian. The advantage of this approach is that the Casimir invariants required for the application of the method are already available, since they only depend on the Poisson bracket. However, this procedure is in general non trivial, because the change to the rotating frame θ ′ = θ -ωt is performed at the level of the single-particle phase space, whereas the Poisson bracket ( 5) is defined on the algebra of functionals of the field variables.

Since the transformation of coordinates to the rotating frame is time-dependent, it is necessary to introduce a new evolution parameter and extend the algebra of observables to the set of functionals U (n e , f, t, k), where t is the time dynamical variable and k its canonically conjugate variable with respect to the extended Poisson bracket

{F, G} t = {F, G} + ∂F ∂t ∂G ∂k - ∂G ∂t ∂F ∂k . (8)
Consequently, we have to consider in the extended phase space the new Hamiltonian H t = H + k. We then look for a new set of dynamical variables (ñ e , f , t, k) such that ñe (r, θ ′ ) = n e (r, θ ′ + ωt),

k = k + W (n e , f ), f (r, θ ′ , v r , v θ ) = f (r, θ ′ + ωt, v r , v θ ), t = t, (9)
where W (n e , f ) is a functional determined by the constraint that the bracket ( 5) is invariant with respect to the change of variables. The mapping (9) determines the following transformation for the derivatives of the functionals:

F χ = Fχ + W χ ∂ F ∂ k , ∂F ∂k = ∂ F ∂ k , ∂F ∂t = ∂ F ∂ t + ω dA ∂ ñe ∂θ ′ Fñe + dAd 2 v ∂ f ∂θ ′ F f , ( 10
)
where χ is either n e or f . Inserting the expressions ( 9)-( 10) into Eq. ( 5), the bracket turns out to be invariant if W = ωP θe +ωP θI , where P θe = -dA(eB 0 /2c)r 2 n e and P θI = dAd 2 v[(q I B 0 /2c)r 2 + M rv θ ]f . As expected, P θe and P θI correspond to the total canonical angular momenta of massless electrons and of the ions, respectively. We restrict the stability analysis to equilibria which are independent on θ. Consequently, from now on we consider physical situations where the gravitational force may be neglected. This approximation also removes the time dependence of the Hamiltonian in the new reference frame. The assumption is valid also for heavy ions but under microgravity conditions [19].

According to the Energy-Casimir method, equilibria of the system correspond to extremals of the free energy functional F given by the linear combination of the Hamiltonian with the Casimirs of the Poisson bracket. Since we consider now a time-independent Hamiltonian in the new reference frame, we omit the contribution of k. For ease of notation we also drop the tilde symbol for the new variables. Stationary equilibria in the rotating frame therefore correspond to extremals of [START_REF] Davidson | An Introduction to the Physics of Nonneutral Plasmas[END_REF]

F = H -ωP θe -ωP θI + C 1 + C 2 , (11)
i.e., to solutions of

F ne = 0, F f = 0. ( 12
)
These are given by

F ′ (n e ) = -h e , G ′ (f ) = -h I , (13)
where the prime denotes the derivative with respect to the argument of the function and where h e = -eφ + ω(eB 0 /2c)r 2 and h I = M v 2 /2 + q I φ -ω[(q I B 0 /2c)r 2 + M rv θ ] indicate the single particle energies in the rotating frame of massless electrons and of the contaminant particles, respectively. Assuming that the arbitrary functions F ′ and G ′ are invertible, the relations [START_REF] Tronci | [END_REF] imply that, at equilibrium, n e = n eeq (h e ) and f = f eq (h I ), where n eeq = (-F ′ ) -1 and f eq = (-G ′ ) -1 . According to the Energy-Casimir method, energy stability (also called formal stability) of the equilibria obtained by setting δF = 0, is guaranteed if δ 2 F , evaluated for such equilibria, has a definite sign for all field variations. We recall that this type of stability implies linear stability (see, e.g. Refs. [17,18])

We consider field variations δn e and δf periodic with respect to θ and vanishing at r = 0 and r = a. We also assume that δf tends to zero for v r and v θ going to infinity. The second variation δ 2 F can then be computed as

δ 2 F (n eeq , f eq ) = dA |∇δφ| 2 8π + F ′′ (n eeq )|δn e | 2 + dAd 2 vG ′′ (f eq )|δf | 2 ,
where ∆δφ = 4π(eδn e -q I d 2 vδf ). The second variation of the free energy functional is then positive if F ′′ (n eeq ) ≥ 0 and G ′′ (f eq ) ≥ 0.

Using Eq. ( 13) and the relations n eeq = (-F ′ ) -1 and f eq = (-G ′ ) -1 , a set of stability conditions is therefore given by dn eeq dh e ≤ 0, df eq dh I ≤ 0. (

In the limit of single-species plasmas, the inequalities (14) separately correspond to the stability condition for a cold pure electron plasma equilibrium in the fluid guidingcenter description and for a kinetic equilibrium of a population of massive ions, respectively [START_REF] Davidson | An Introduction to the Physics of Nonneutral Plasmas[END_REF]20,21]. Defining ne (r) = n eeq (h e ), f (r, v r , v θ ) = f eq (h I ) and nI (r) = d 2 v f (r, v r , v θ ). Conditions ( 14) imply for confined equilibria with monotonically decreasing density profiles (n ′ e (r) ≤ 0, n′

I (r) ≤ 0) that ω ≥ 4π r 2 c B 0 r 0 dss(en e (s) -q I nI (s)) , (15)
and

ω -≤ ω ≤ ω + , (16)
respectively, with

ω ± = ω ci 2 1 ± 1 + 16πq I M ω 2 ci r 2 r 0 dss(en e (s) -q I nI (s)) 1/2 , (17)
where ω ci = |q I |B 0 /M c is the ion cyclotron frequency and use has been made of the relation φ ′ (r) = (4π/r) r 0 ds(en e (s) -q I nI (s)), obtained from Eq. (3) imposing the condition φ ′ (0) = 0.

Close to the axis (so that ne (r) ≈ ne (0) and nI (r) ≈ nI (0)), the reality condition for ω ± corresponds to the self-field parameter 4(ω D I + ω D e )/ω ci to be less than unity, with ω D I = 2π|q I |cn I (0)/B 0 and ω D e = 2πecn e (0)/B 0 , the ion and electron diocotron frequencies, respectively. In particular, if 4(ω D I +ω D e )/ω ci ≪ 1, the conditions (15) and ( 16) reduce, at the leading order, to

ω ≥ ω D e + ω D I , (18)
ω D e + ω D I ≤ ω ≤ ω ci . (19)
Clearly, the condition ( 18) is implied by the condition (19). The latter relation identifies a finite range of frequencies ω for which the stability condition can be satisfied for confined equilibria. When the ion species acts as a small contamination of a pure electron plasma, the condition 4(ω D I + ω D e )/ω ci ≪ 1 reduces to 4ω D e ≪ ω ci , i.e., ne (0) ≪ (1/8π)|q I |B 2 0 /eM c 2 . This corresponds to a modified Brillouin electron density limit, which depends also on the charge and the mass of the contaminant species. In the opposite case, ne (0) = 0, the usual Brillouin density limit for the negative ion species is retrieved, nI (0) ≪ (1/8π)B 2 0 /M c 2 . An important example of stable confined equilibrium can be obtained from an electron density and an ion distribution function corresponding to n eeq (h e ) = ne (0) exp(-h e /κT ) and f eq (h I ) = (n I (0)M/2πκT ) exp(-h I /κT ), respectively, where κ denotes the Boltzmann constant and T is a constant temperature. For this case, the conditions ( 14) are clearly satisfied and the confinement condition (19) follows requiring n′ I (r) = (n I (0)/κT )(-ω 2 M r + ω|q I |B 0 r/c -|q I |φ ′ ) exp(-h I /κT ) to be negative. The stability and confinement of non-neutral plasmas in thermal equilibrium have been discussed in Ref. [22] and confinement conditions analogous to Eq. ( 19) for kinetic multi-species thermal plasmas have been described in Ref. [START_REF] Davidson | An Introduction to the Physics of Nonneutral Plasmas[END_REF]. In the hybrid model presented here, where the electron species is treated as a massless fluid, for stable equilibria with small self-field parameter the confinement condition for the negative ion species is sufficient to guarantee the con-finement of both species.

This work was supported by the European Community under the contract of Association between EURATOM, CEA, and the French Research Federation for fusion studies. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Financial support was also received from the Agence Nationale de la Recherche (ANR GYPSI) and from the CNRS PEPS project GEOPLASMA. This work was also partially supported by the Italian Ministry for University and Research (MIUR) "PRIN 2009" funds (Grant No. 20092YP7EY). E. T. acknowledges fruitful discussions with Eric Soccorsi as well as useful interactions with Cesare Tronci and Phil Morrison about Hamiltonian hybrid models.


E Tassi

D Grasso

L Comisso

Linear stability analysis of collisionless reconnection in the presence of an equilibrium flow aligned with the guide field

The influence of a velocity jet, directed along a magnetic guide field, on the linear evolution of collisionless reconnection is investigated both analytically and numerically. The analysis covers both the small and large ∆ ′ regimes, with ∆ ′ indicating the standard tearing stability parameter, and is carried out, in slab geometry, by means of a reduced four-field model for magnetic reconnection accounting for two-fluid effects. Analytical dispersion relations are derived in both regimes and their predictions on the growth rates are tested against numerical simulations. In both regimes the presence of the flow is shown to have a stabilizing effect, with growth rates decreasing when increasing the amplitude of the equilibrium flow. The analytical results predict that a decrease in the growth rate could be obtained also by reducing the characteristic width of the equilibrium flow profile. Such stabilizing effects appear to be stronger in the small ∆ ′ regime. A very good quantitative agreement is found between the analytical predictions and the numerical results. As a complement to the analysis, we also consider, in the small ∆ ′ regime, the dispersion relation in the absence of equilibrium flow, which extends a previously derived dispersion relation by including a corrective term due to plasma parallel compressibility. It is shown that such correction can have a stabilizing effect and yields a better agreement with the numerical results.

I. INTRODUCTION

Understanding the interplay between magnetic reconnection events and shear flows is of great relevance for plasmas both in fusion and astrophysical contexts. In tokamak plasmas velocity flows are, for instance, generated by the neutral beam injection [1] or ion cyclotron heating [2] and their coexistence with tearing modes is of great importance for the experimental performance. On the other hand, in space plasmas large scale jets are often observed in regions where a magnetic field shear is also present, such as at the Earth's magnetopause or in the solar corona [3].

In investigating the effects of flows on reconnection events, one of the immediate and natural questions concerns the impact that such flows can have on the linear stability and on the growth rate of reconnective perturbations. Indeed, early analytical studies on this subject date back already to some decades ago. In particular, most of these linear stability studies were concerned about the role of shear flows, in general lying in the same plane of the sheared equilibrium magnetic field. Effects of such flows on the resistive tearing instability were investigated studying the interplay of the latter with the Kelvin-Helmholtz instability [4,5]. The dependence of the linear stability threshold and growth rate on the intensity of the planar shear flow was studied analytically and numerically in Refs. [6][7][8][9]. On the other hand, the stabilizing effect of a linear flow on a current sheet in two-dimensions was discovered in Ref. [10]. In cylindrical geometry, the influence of an equilibrium flow on the tearing stability parameter ∆ ′ was investigated in Ref. [11] and, recently, accounting for more general flows and finite pressure effects, in Ref. [12].

In the present paper we consider, in slab geometry, the influence of an out-of-plane flow (i.e. directed perpendicular to the plane of the sheared equilibrium magnetic field) on the linear phase of collisionless reconnection. We remark that, recently the role played by an out-of-plane velocity jet on the fast reconnection process in magnetohydrodynamics (MHD) and Hall-MHD has been highlighted and investigated by means of nonlinear simulations, showing that it can significantly modify the structure of the out-of-plane magnetic field and change the rate of the reconnection process [13,14]. Also, in the tokamak context, the stabilizing role of a purely toroidal flow on classical and neoclassical tearing modes has been highlighted [16] and its influence on ∆ ′ has been determined analytically [17]. Clearly, the latter two works account also for toroidal effects that are absent in the slab geometry that we consider in the present paper.

Our stability analysis is based on a four-field model [18] where reconnection is mediated by electron inertia and the presence of a strong guide field is assumed. We point out that the above mentioned works on linear stability, either did not consider out-of-plane velocity equilibria [10,11] or they included them as special cases and adopted three-dimensional (3D) incompressible MHD as model equations [4][5][6][7][8]12]. The model we consider here, on the other hand, is reduced to two dimensions but it retains two-fluid and compressibility effects. In the presence of such effects, even in the 2D approximation, the presence of an equilibrium out-of-plane flow can in principle modify the linear evolution of magnetic reconnection (in pure incompressible 2D MHD, on the contrary, an out-of-plane equilibrium velocity field such as the one we consider here, would have no impact on linear stability). Also, whereas previously mentioned works (with the exception of Ref. [9]) focused on dissipative resistive tearing modes, here we consider inertial reconnection, which is an intrinsically non-dissipative process. The adopted model has indeed been shown to correctly possess a Hamiltonian structure [19].

In order to carry out the linear stability analysis, we make use of standard asymptotic matching techniques and, following previous linear works [21][22][23], we derive a dispersion relation providing the growth rate of the collisionless tearing instability in the presence of an out-of-plane velocity jet, symmetric with respect to the resonant surface. As paradigmatic case we consider a velocity profile corresponding to a Bickley jet, which can model intense flows and which has been previously adopted for instance to investigate the interplay between the Kelvin-Helmholtz instability and the tearing instability in the reconnection plane [24].

We derive a dispersion relation in both the small and large ∆ ′ regimes. Subsequently, in order to validate the analytical results, we check them against numerical simulations. In particular, for the quantitative analysis, we consider flow amplitudes from the sub-Alfvénic range, relevant for laboratory plasmas, to the Alfvénic range, occurring in space plasmas.

We also consider different values of the plasma β, where β = (5/3)µ 0 p 0 /B 2 0 is proportional to the ratio between the background kinetic pressure and the toroidal magnetic pressure (we adopted the definition of β given in Ref. [18], which explains the unusual factor 5/3 coming from the choice of an adiabatic equation of state) . Such parameter is indeed crucial for the effect of the out-of-plane flow on the dynamics.

Our analysis is primarily meant to help understanding basic properties of collisionless reconnection. Given the relative simplicity and generality of the adopted model, we do not claim that we can make very detailed predictions on a specific experimental situation or observation. However, our choice for the orderings and for most of the parameter values in the following, will be suggested by typical situations occurring in tokamaks which, as above mentioned, provide one of the main motivations for studying the role of out-of-plane flows on reconnection. Moreover, it is namely in weakly collisional plasmas, such as tokamak plasmas, that reconnection mediated by electron inertia is believed to play an important role [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF]. We will consequently focus mainly on the low-β regime (i.e. β ≪ 1), although the model is in principle valid also for higher β values. Evident limitations of the model, with respect to an application to tokamaks, come from the use of a slab geometry and the absence of kinetic effects. On the other hand, such reduced model is amenable to a fully analytical treatment, and might provide a leading-order indication about the role of a class of out-of-plane (i.e. "toroidal") flows in collisionless reconnection with two-fluid effects. Our choice for the equilibrium velocities also presents some simplifications with respect to the generic situation in a tokamak. In particular it assumes the absence of shear and a nonmonotonic velocity profile around the resonant surface. Non-monotonic velocity profiles, on the other hand, have been adopted in order to show the stabilizing effect of toroidal flows against sawteeth in tokamaks [START_REF] Kleva | [END_REF]27]. Also, even in the absence of shear, toroidal flows, have been shown to have a stabilizing effect [16]. Thus, we consider the present analysis as the investigation of the role of a particular class of out-of-plane flows, which might be relevant, but that will need to be generalized in the future, in particular with the inclusion of shear on the resonant surface and asymmetries.

We anticipate that, on the basis of the results that we obtained, an out-of-plane shear flow with the same direction of the guide field turns out to have a stabilizing effect on the tearing mode, whose growth rate decreases with increasing amplitude of the equilibrium flow. This stabilizing effect is more evident in the small ∆ ′ regime and adds to the one following from the plasma compressibility in the toroidal direction and that was neglected in a previous analysis carried out in absence of shear flow [18].

The paper is organized as follows: section II describes the model and its linearized version; in sections III and IV the dispersion relations valid in the small and large ∆ ′ regimes respectively are derived; in both these two sections the comparison with the results obtained from numerical simulations are also presented. Conclusions are drawn in the last section.

II. REVIEW OF THE MODEL EQUATIONS AND LINEARIZATION

We consider the following two-dimensional model for non-dissipative reconnection mediated by electron inertia: [18]

∂(ψ -d 2 e ∇ 2 ψ) ∂t + [φ, ψ -d 2 e ∇ 2 ψ] -d β [ψ, Z] = 0, (1)
∂∇ 2 φ ∂t + [φ, ∇ 2 φ] + [∇ 2 ψ, ψ] = 0, (2)
∂Z ∂t + [φ, Z] -c β [v, ψ] -d β [∇ 2 ψ, ψ] = 0, (3)
∂v ∂t + [φ, v] -c β [Z, ψ] = 0. (4)
Eqs. ( 1)-(4) determine the evolution of the magnetic field 4) is written according to the following normalization:

B = ∇ψ × ẑ + (c β Z + 1)ẑ,
t = t t A , x = x L , ψ = ψ LB 0 , Z = Ẑ B 0 , φ = φ v A LB 0 , v = v v A , (5)
where carets denote dimensional quantities. In (5) B 0 indicates the amplitude of the magnetic guide field, L is the equilibrium magnetic field characteristic scale length, whereas t A and v A are the characteristic Alfvén time and speed, based on B 0 . In the model equations, Eq. ( 1) is an Ohm's law for a collisionless plasma, accounting for electron inertia. Eq. ( 2) is a vorticity equation. Eq. ( 3) governs the evolution of the out-of-plane perturbation of the magnetic field, but also the electron pressure, given that the derivation of the model assumes p e ∝ -c β Z, with p e indicating namely the electron pressure perturbations. Eq. ( 4) is an evolution equation for the out-of-plane component of the velocity of the plasma center of mass, which corresponds approximately to the ion fluid velocity. We recall that the derivation of the model also assumes that an approximate magnetohydrodynamic equilibrium is maintained in the perpendicular plane. This corresponds to filter our the compressional Alfvén waves and leads to the above proportionality relation between pressure and magnetic field perturbations. An adiabatic equation of state for the electron pressure is assumed, which permits, together with the evolution of equation for the out-of-plane magnetic perturbations and with the above mentioned proportionality relation, to close the system by eliminating the divergence of the electron velocity field, which is assumed to be small but finite. Such compressibility effect is retained thanks to the assumption of finite β and eventually leads to the appearance of the coefficient c β . In the limit case β = 0, Eqs. ( 3)-( 4) decouple from the system and the two remaining equations correspond to the equations of reduced magnetohydrodynamics with electron inertia corrections.

The model is derived adopting a drift approximation, which assumes frequencies lower than the ion cyclotron frequency and a strong guide field along the z direction. As above mentioned, the ion and electron fluids are assumed to be incompressible at the leading order, although compressibility effects are retained in the evolution equations for the electron pressure and the out-of-plane perturbations of the magnetic field. No ion pressure effects are taken into account, on the other hand, for the ion fluid is assumed to be cold. The derivation of the model can be found in Ref. [18].

Consistently with the conservative two-fluid system, from which it is derived, the system (1)-( 4) possesses time-reversal symmetry with respect to the transformation ( t, x, B, v, φ) → (-t, x, -B, -v, φ). Also, as shown in Ref. [19], the model possesses four infinite families of Casimir invariants, reflecting, among other properties, the conservation of the parallel canonical momentum of the ion fluid, and of the generalized vorticity.

We consider the system on the domain {(x, y) : -∞ ≤ x ≤ +∞, -L y ≤ y ≤ L y }, with constant L y , and linearize it about the following equilibrium solutions:

ψ = ψ eq (x), φ = 0, Z = 0, v = v eq (x), (6)
with the functions ψ eq (x) and v eq (x) to be specified. The functions indicated in Eq. ( 6)

are evidently equilibrium solutions of the four-field model ( 1)-( 4). However, this choice for the equilibria clearly excludes important phenomena, such as for instance drift waves, which are associated with the presence of equilibrium pressure gradients, or fluid-like instabilities due to the presence of an equilibrium E × B velocity. Given that the focus of the present paper is on the effect of flows directed along the guide field, we carry out our analysis in this simplified frame and intend to study more general equilibria in a future work.

We require the planar component of equilibrium magnetic field B eq (x) = -dψ eq /dx to vanish at x = 0, the latter corresponding then to a resonant surface.

We consider perturbations of the form f (x, y, t) = f (x) exp(γt + iky) + c.c, where f is a generic field perturbation, whereas γ and k indicate the corresponding frequency and wavenumber.

The linearized system then reads:

γ[-d 2 e ψ′′ + (1 + d 2 e k 2 ) ψ] -ik(ψ ′ eq -d 2 e ψ ′′′ eq ) φ -ikd β ψ ′ eq Z = 0, (7)
γ( φ′′ -k 2 φ) + ikψ ′′′ eq ψ -ikψ ′ eq ( ψ′′ -k 2 ψ) = 0, (8)
γ Z -c β ik(v ′ eq ψ -ψ ′ eq ṽ) -d β ik[ψ ′′′ eq ψ -ψ ′ eq ( ψ′′ -k 2 ψ)] = 0, (9)
γṽ -ikv ′ eq φ + c β ikψ ′ eq Z = 0, ( 10
)
where the prime denotes the derivative with respect to x. We assume, as in standard linear theory for tearing modes, the following symmetries

ψ(x, -y) = ψ(x, y), v(x, -y) = v(x, y), (11)
φ(x, -y) = -φ(x, y), Z(x, -y) = -Z(x, y), (12)
which are respected by the model equations ( 1)-( 4). As a consequence, in ( 7)-( 10), ψ and ṽ are purely real-valued functions, φ and Z are purely imaginary, and the growth rate γ is purely real.

Note that the effect of the equilibrium velocity on the fields ψ and φ enters through the term depending on d β in Eq. (7). Remark that d β reflects a two-fluid nature of the model.

In particular it is worth recalling that d β reduces to ρ s , i.e. the Larmor radius of ions with electron temperature, for c β ∼ √ β.

We assume now perturbations vanishing at x → ±∞ and consider the following choice for the equilibrium solutions:

ψ eq (x) = -a ln cosh x, v eq (x) = v cosh 2 x , (13)
with constant a and v. We anticipate, however, that it is straightforward to extend our analysis to other standard choices of magnetic equilibria for the tearing stability problem (e.g. ψ eq (x) = a cos x or ψ eq (x) = a/ cosh 2 x and analogously for the velocity equilibria), provided that some symmetries of the system (1)-( 4) are respected. In particular, we look for even solutions for ψ and v and odd solutions for φ and Z.

The equilibrium (13) describes a sheared magnetic field with a resonant surface at x = 0 and an out-of-plane velocity jet.

As is customary for collisionless tearing stability analysis, we consider modes that evolve slowly on the characteristic Alfvén time and a reconnection layer thickness (which is typically of the order of d e ) much smaller than the equilibrium scale length:

γ ≪ 1, d 2 e ≪ 1. (14)
The existence of a resonant surface requires, however, the boundary problem ( 7)-( 10) to be solved separately in a thin inner layer surrounding the resonant surface, where the terms proportional to d 2 e cannot be neglected, and in one, or more, outer regions. The matching condition between the inner and outer solutions will provide the relation determining the growth rate of the tearing mode as function of the parameters of the system.

In our analysis we consider two regimes, depending on the value of the ∆ ′ parameter, which is defined as

∆ ′ = lim x→0 + 1 ψ d ψ dx -lim x→0 - 1 ψ d ψ dx , (15)
where ψ refers here to the solution in the outermost region.

In the following sections we derive, under the above hypotheses, dispersion relations valid for the equilibrium (13), treating separately the cases ∆ ′ d e ≪ 1 (small ∆ ′ regime) and ∆ ′ d e ≫ 1 (large ∆ ′ regime).

III. STABILITY ANALYSIS: SMALL ∆ ′ REGIME

In the small ∆ ′ regime the analysis can be carried out by introducing an inner region, with a size of the order of d e , centered around the resonant surface, and an outer region, corresponding to the complement of the inner region.

A. Outer region

Neglecting the terms proportional to the electron inertia, we obtain from ( 7)-( 10) that in the outer region the relevant system is given by

γ ψ -ikψ ′ eq φ -ikd β ψ ′ eq Z = 0, (16)
γ( φ′′ -k 2 φ) + ikψ ′′′ eq ψ -ikψ ′ eq ( ψ′′ -k 2 ψ) = 0, ( 17
) γ Z -c β ik(v ′ eq ψ -ψ ′ eq ṽ) -d β ik[ψ ′′′ eq ψ -ψ ′ eq ( ψ′′ -k 2 ψ)] = 0, ( 18
)
γṽ -ikv ′ eq φ + c β ikψ ′ eq Z = 0. ( 19
)
Neglecting the terms proportional to γ in (17), we obtain, for the equilibrium (13) [20]:

ψ(x) =      ψ 0 e -kx 1 + tanh x k , if x > 0 ψ 0 e kx 1 -tanh x k , if x < 0 (20)
where ψ 0 is a constant. We remark that, from such solutions, one obtains

∆ ′ = 2 1 k -k . ( 21
)
Then we assume Z ∼ 0 in the outer region. This assumption leads us to the reduced MHD ordering, for the perturbations of the magnetic flux function and of the stream function in the outer region. This simplification allows to obtain a complete solution for the outer system providing the correct matching with the inner solution, while retaining the effects of the equilibrium velocity field in the outer region, which we do not want to bound a priori.

Also, we do not impose a priori restrictions on the value of d β .

From ( 16) then we obtain

φ ∼ -i γ kψ ′ eq ψ, (22)
and finally, from (22) and ( 19) one finds

ṽ ∼ v ′ eq ψ ′ eq ψ. ( 23
)
We observe that in the vicinity of the resonant surface, the behavior of the outer solution for ψ is given by

ψ = ψ 0 1 + 1 2 ∆ ′ |x| + O(x 2 ) . (24)
In the small ∆ ′ regime that we are considering, the outer solution can then be considered as a constant when approaching the resonant surface.

B. Inner region

In the inner region, the terms proportional to the electron inertia are retained, but the following local expansions

ψ ′ eq ∼ -ax, ψ ′′′ eq ∼ 2ax, v ′ eq ∼ -2vx, (25)
are adopted for the equilibria. The resulting system reads

-d 2 e γ ψ′′ + γ(1 + d 2 e k 2 ) ψ + ikax(1 + 2d 2 e ) φ + ikd β ax Z = 0, (26)
γ φ′′ -γk 2 φ + ikax[ ψ′′ + (2 -k 2 ) ψ] = 0, ( 27
) γ Z -ikd β ax ψ′′ + 2 -k 2 -2 v d i a ψ -ikc β axṽ = 0, ( 28
)
γṽ + ik2vx φ -ikc β ax Z = 0. ( 29
)
In the inner layer the terms involving derivatives along the x direction dominate over those with derivatives along y. Therefore we have k 2 ψ ≪ ψ′′ and k 2 φ ≪ φ′′ . More in general, only the leading order x-derivative terms should be retained in the inner region [22]. This would imply, in general, the suppression of all the terms proportional to ψ in Eq. ( 28). In particular, in the presence of a moderate equilibrium flow, with v ∼ ad i , the contributions due to the equilibrium velocity would be negligible in the inner region and already available dispersion relations [18,22] might apply. For stronger out-of-plane flows with v ≫ ad i , however, such terms can in principle affect the growth rate of the perturbation. Moreover, we observe that, if we generalize our velocity equilibrium by taking v eq (x) = v/ cosh 2 (x/l) with l arbitrary positive constant, the outer solutions remain unchanged and also the analysis of inner layer is the same, upon replacing v with v/l 2 . Therefore, the contributions of the equilibrium flow do not necessarily require a very large velocity to be relevant, but they can influence the inner layer equations also if the jet is sufficiently narrow, so that v ≫ l 2 ad i . In the inner layer we keep only the leading order x-derivative terms in ( 26)-( 29) but, considering the above discussed regimes, we retain also the term due to the equilibrium flow in Eq. ( 28), although it does not involve x-derivatives.

The resulting system can be written in the following way:

gd 2 e ψ′′ -g ψ -ix φ -id β x Z = 0, (30)
g φ′′ + ix ψ′′ = 0, ( 31
) g Z -id β x ψ′′ + ic β υx ψ -ic β xṽ = 0, ( 32
)
gṽ + iυx φ -ic β x Z = 0, (33)
where

g = γ ak , υ = 2 v a . ( 34
)
Rearranging the terms and neglecting d 2 e when compared to d 2 i , one obtains from ( 30)-(33), the following relation:

ψ + i x g φ = g 2 d 2 e + d 2 β x 2 g 2 + αc 2 β x 2 ψ′′ , (35)
where α = 1 + υd i . We observe that, in terms of dimensional quantities, one in general

has υd i = 2( v/ω cipol L)
, where v is the dimensional amplitude of the velocity equilibrium, whereas ω cipol indicates the ion cyclotron frequency based on the characteristic equilibrium "poloidal" magnetic field B pol = aB 0 .

At this point one can take advantage of the small ∆ ′ hypothesis, and assume that ψ is approximately constant around the resonant surface [28]. More precisely, we set ψ = ψ0 + ψ1 (x), with | ψ1 | ≪ ψ0 . If we choose the value of the arbitrary constant ψ0 equal to one, and neglect ψ1 when compared to 1, from Eq. ( 35) we obtain

ψ′′ 1 = g 2 + αc 2 β x 2 g 2 d 2 e + d 2 β x 2 (1 -x ξ), (36)
where we introduced the displacement function ξ = -i φ/g.

The matching between the outer and inner solution can be performed at the level of the ∆ ′ parameter. Taking advantage from the constant-ψ approximation, one can make use of the relation

∆ ′ = +∞ -∞ dx ψ′′ , (37)
which connects the limit for large x of the inner solution ψ(x) with ∆ ′ , the latter carrying information about the outer solution. Combining (36) and (37), we find that the matching condition can be written as

∆ ′ = +∞ -∞ dx g 2 + αc 2 β x 2 g 2 d 2 e + d 2 β x 2 (1 -x ξ). ( 38
)
Determining the inner solution for the displacement function will provide then, through Eq.

(38) a relation between the growth rate g and the various parameters of the system, among which ∆ ′ , which is known from the outer solution, according to (15).

We observe that the effect of an out-of-plane equilibrium flow, whose signature is the α parameter, combines in Eq. ( 38) with the c 2 β contribution, indicating that the presence of the flow can modulate the effects due to the plasma β. On the other hand, in the limit c β → 0, the out-of-plane flow dynamics decouples from the system, which is reflected by the fact that in that limit, no influence of the flow would be present in the relation (38).

Analogously to Ref. [22], we introduce at this point the parameters

δ = gd β c β 1/2 , σ = g c β d β , ǫ = d 2 e d 2 i σ, (39)
and the following renormalized dependent and independent variables:

ξ = δ ξ, z = x δ . ( 40
)
With the help of ( 31) and (36), we obtain that the renormalized displacement function has to satisfy the equation

ξ′′ -z αz 2 + σ z 2 + ǫ (z ξ -1) = 0. ( 41
)
In the limit of negligible equilibrium flow (which corresponds to α = 1), Eq. ( 41) reduces to the corresponding inner layer equation derived in Ref. [22]. Requiring the antisymmetry of φ and a convergence toward the outer region, leads to the following boundary conditions for (41):

ξ(0) = 0, lim z→±∞ ξ(z) = 0. (42)
Because of the symmetry of the problem, however, we consider only z ≥ 0.

We assume the ordering

ǫ ≪ σ α ≪ 1. ( 43
)
Given that d 2 i /d 2 e ≫ 1, the inequality ǫ ≪ σ/α can be approximated and rewritten as υ ≪ d i /d 2 e , or in dimensional terms, as v ≪ (M i /m e )Lω cipol , where M i and m e indicate the ion and electron mass, respectively. This ordering can apply to various laboratory plasmas and in particular to tokamak plasmas . We anticipate that our analysis holds in principle also for the alternative ordering ǫ ∼ σ/α. However, very narrow velocity profiles v eq (x) = v/ cosh 2 (x/l), with l such that v ∼ l 2 ad i /(2d 2 e ), and/or super-Alfvénic equilibrium flows, would be required in order to reach ǫ ∼ σ/α. These situations are unlikely to occur in the laboratory plasmas.

The inequality σ/α ≪ 1, or, equivalently, γ ≪ akc β d β (1+υd i ), on the other hand, applies to finite β plasmas near marginal stability.

The ordering (43) allows us to solve Eq. (41) in two separate regions. For z → 0, the leading order solution of Eq. ( 41) is ξ(z) ∼ 0.

(44)

On the other hand, for z 2 ≫ σ/α, Eq. ( 41) can be approximated by ξ′′ -αz 2 ξ = -αz.

(45)

This boundary value problem can be mapped into the one solved in Ref. [21], by rescaling the dependent and independent variables by -α -1/4 and α 1/4 , respectively. The solution at large z for our case is then given by

ξ(z) ∼ α 1/2 2 z 1 0 dt(1 -t 2 ) -1/4 exp -α 1/2 tz 2 /2 . ( 46
)
With the help of the approximate solutions for ξ, we can now evaluate the expression (38), which we rewrite as

∆ ′ = δc 2 β d 2 β αI 1 + σ -ǫ ǫ -(α -1) I 2 , (47)
where

I 1 = +∞ -∞ dz(1 -z ξ), I 2 = +∞ -∞ dz 1 -z ξ 1 + z 2 /ǫ . ( 48
)
Because σ/α ≪ 1, the dominant contribution in I 1 turns out to be that given by the solution (46) for z 2 ≫ σ/α. Therefore we get [29]

I 1 ≃ 2π α 1/4 Γ(3/4) Γ(1/4) ≃ 2.12 α 1/4 , ( 49
)
where Γ indicates the Gamma function. On the other hand, considering the ordering (43), one can show that the dominant contribution in I 2 is given by

I 2 ≃ 2ǫ σ/α 0 dz z 2 + ǫ ≃ πǫ 1/2 . ( 50
)
Combining the results we finally obtain the following dispersion relation

d β ∆ ′ = α 3/4 c 3/2 β d 1/2 β 2.12 + g 1/2 π d e 1 -αd 2 e c 2 β d 2 β g 1/2 . ( 51
)
For positive ∆ ′ , Eq. ( 51) always possesses one positive solution for g 1/2 . In particular we remark that the growth rate vanishes for ∆ ′ = 0, as in standard tearing instability.

In the absence of equilibrium flow, Eq. ( 51) reduces, up to corrections of order d 2 e /d 2 i , to the dispersion relation of Ref. [22] without the terms due to diamagnetic drift and resistivity.

The effects due to the equilibrium flow enter the dispersion relation in two places, and in both cases, not surprisingly, they get suppressed if c β = 0. The coefficient α 3/4 tends to amplify the effect of a finite plasma β, for an equilibrium bulk velocity flowing in the direction of the guide field. The second contribution of the plasma flows, on the other hand, couples also with the electron skin depth, which is related to the mechanism breaking the frozen-in condition in this model. We remark, however, that for finite c β , the relation

1 -αd 2 e c 2 β /d 2 β ≈ 1 -2(v/a)(d 2 e /d i ) holds.
Therefore, if we compare the magnitude of the last two terms in (51), we obtain that the latter modification due to the equilibrium flow is generally small, unless extremely high velocities (v ∼ ad i /(2d 2 e ) ) are attained, or very narrow jets (l ∼ (2d 2 e v/ad i ) 1/2 ) are considered. Actually, in order to be consistent with the first inequality in (43), such modification due to the equilibrium flow should be neglected. However, as above anticipated, the dispersion relation (51) applies in principle also in the

regime α ∼ d 2 i /d 2 e .
It is also worth noticing that, in general, the tendency of a field-aligned equilibrium flow, appears to be that of reducing the effect of the term depending on the electron skin depth in the dispersion relation. Incidentally, we remark that the condition v ∼ ad i /(2d 2 e ) is namely the condition that, in Eq. ( 28), makes the term related to the equilibrium velocity of the same size as the term proportional to ψ′′ , assuming that in the inner layer d 2 /dx 2 ∼ 1/d 2 e . In order to test the validity of the dispersion relation (51), we compared the predictions of the growth rate obtained from (51), with the results of numerical solutions of the four-field model. The numerical code we adopted for the simulations solves the four-field system with initial conditions obtained by perturbing the equilibria (6), with ψ eq (x) and v eq (x), given by (13). The model equations are solved on a grid consisting of up to 1024×512 grid points.

All the fields are split in the time independent equilibrium and an evolving perturbation, which is advanced in time according to a third order Adams-Bashforth algorithm. Periodic boundary conditions have been imposed along the y direction, whereas Dirichlet conditions have been applied along the x direction imposing that all the perturbed fields vanish at the boundaries. A pseudospectral method is adopted for the periodic direction, while a compact finite difference algorithm on a non-equispaced grid is used for the spatial operations along the x direction. The tearing instability is initiated by perturbing the equilibrium with a small disturbance on the parallel current density j = -∇ 2 ψ of the form δj (x, y) = δj (x) cos(2πy/L y ), where δj (x) is a function localized within a width of the order d e around the rational surface x = 0. The numerical growth rates have been determined by evaluating, during the linear phase, the value of d/dt log |δψ X (t)|, with δψ X (t) indicating the value of ψ at the X-point of the magnetic island. For the linearized system this quantity corresponds to the growth rate of the m = 1 mode. Fig. 1 shows a comparison between analytical and numerical results, over a wide range of values of flow amplitude υ and for two values of c β . We observe that the growth rates decrease when increasing υ, indicating that the out-of-plane equilibrium flow has a stabilizing effect on the reconnective perturbation. The analytical growth rates agree quantitatively well with those obtained from the numerical simulations, reproducing the same tendency. The relative discrepancies between analytical and numerical values are typically well below 10 %, in particular for sub-Alfvénic or slightly super-Alfvénic flows, which are those of interest for most applications. We observe also that the tendency toward stabilization shown by the numerical results, is caught by our dispersion relation, even for v ∼ ad i , where, in principle, contributions from ψ ′′′ eq (x), which we discarded in the inner region, might have played a role. It is also opportune to compare the dispersion relation (51), in the limit α = 1, with the dispersion relation which was derived in Ref. [18] for the same four-field model but with a purely magnetic equilibrium in the small ∆ ′ regime. The latter relation differs from ours because it does not account for the first term on the right-hand side of (51), which is due to the contribution of I 1 . Such contribution originates from plasma compressibility in the direction parallel to the guide field.

FIG. 2. Comparison between growth rates obtained from the dispersion relation of Ref. [18] (blue dash-dotted line), from the dispersion relation (51) (red dashed line) and from numerical simulations (diamonds connected by a solid line), for different values of c β . The figure indicates that the dispersion relation (51) yields a better agreement with respect to the dispersion relation of Ref. [18], due to the presence of a stabilizing corrective term related to the plasma parallel compressibility. The values of the parameters are ∆ ′ = 0.9, d e = 0.2, d i = 8.6, a = 1, v = 0 and the simulation domain is the same as in Fig. 1.

In Fig. 2 we compare growth rates obtained from numerical simulations with those determined with our relation (51) and with the dispersion relation of Ref. [18]. We observe that the correction due to I 1 yields a better agreement with the numerical results. The correction appears to become relevant for large values of c β , corresponding to β 1. In particular, this suggests that this correction is unlikely to be relevant for tokamak plasmas. Nevertheless we remark that Fig. 2 shows that the corrective term yields slightly smaller growth rates, thus indicating a stabilizing effect due to the finite plasma parallel compressibility.

where we made use of Eq. (33) to eliminate ṽ from the system. From (54)-(56) we can obtain the relation

Z = - d β 1 + α c 2 β g 2 x 2 φ′′ . (57)
Differentiating (54) twice with respect to x and combining it with (55) and (57) leads to the following equation

  x 2   d 2 β 1 + α c 2 β g 2 x 2 φ′′   ′ -(g 2 + x 2 ) φ′   ′ = 0. ( 58
)
The solution of (58) for large x, has to match the ideal MHD solution (53). Adopting the procedure of Ref. [18] we introduce

s = x g , V = g φ 0 1 1 + α c 2 β g 2 x 2 d 2 φ ds 2 , ( 59
)
where the function φ(x), when transformed into a function of s, has been renamed with the same symbol. Making use of the large ∆ ′ hypothesis, Eq. ( 58), integrated once, can then be reduced to d ds

s 2 1 + s 2 dV ds -Q 2 (1 + αc 2 β s 2 )V = 0, ( 60
)
where Q = g/d β . In the limit s → 0, the leading order solution, which possesses the desired behavior at large x, is given by

V = V 0 [s -1-ν + A i s ν + O(s 1-ν )], (61)
where V 0 is an arbitrary constant and

ν = 1 2 1 + 4Q 2 -1 , A i = 1 Q G Q α 1/2 c β . (62)
In (62) the function G(x) is defined by

G(x) = √ x 2 Γ(1/4 + x/4) Γ(3/4 + x/4) . ( 63
)
The expression for ν in (62) comes as a solution of the indicial equation that one obtains when solving (60), in the limit s → 0, in terms of power-series solutions. On the other hand, the constant A i gets determined by matching the solution for s → 0 with a solution of (60), valid for s ≫ 1, through an intermediate matching (details for this procedure can be found in Ref. [18]).

C. Electron region

In the innermost region, both terms proportional to d e and d β are retained. However, we assume that the shear Alfvén wave, the sound wave and the equilibrium flow play no role in this region, which is mainly governed by electron dynamics. More precisely, we assume

g ≫ c β d e , d 2
e d i υ ≪ 1, g ≫ d e d 2 e d i υ 1/2 , (64)
which allows us to neglect the above mentioned effects. Indeed, from Eqs. ( 32)-(33), one obtains the relation

(g 2 + c 2 β x 2 ) Z = id β gx ψ′′ -ic β gυx ψ + c β υx 2 φ. ( 65
)
Given that we are considering the electron region, we can make the approximations x ∼ d e and d/dx ∼ 1/d e . Neglecting the effect of sound waves amounts to considering c 2 β x 2 ≪ g 2 , which yields the first inequality in (64). The contributions coming from the equilibrium velocity can be discarded in the electron region, if the second and third term on the righthand side of Eq. ( 65) are negligible with respect to the first term on the right-hand side. This yields the second and third inequalities in (64), respectively, where, for the latter inequality we also made use of Eq. (31). We remark that the second inequality in (64) essentially corresponds to the inequality ǫ ≪ σ/α assumed in (43) for the small ∆ ′ case. Once that Eq. ( 65) is simplified by means of the above assumptions, the term -ix φ in (30), which is the one related to shear Alfvén waves, turns out to be negligible due to the assumption d β ≫ d e .

It is convenient, at this point, to apply the following Fourier transform along the

x direction f (p) = +∞ -∞ dx f (x)e -ipx , (66)
to a generic perturbation f (x). The system is then mapped into

g(1 + d 2 e p 2 ) ψ = d 2 β d dp (p 2 φ), ( 67
)
gp 2 φ = d dp (p 2 ψ), (68)
Z = d β g d dp (p 2 ψ). (69)
Upon introducing the rescaled quantity r = d e p, from (67)-(69), one obtains the equation admits g = 0 as solution, again indicating complete stabilization but only in an asymptotic limit.

d dr r 2 1 + r 2 d Z dr -Q 2 Z = 0, (70)
The quantitative agreement with the numerical predictions is again satisfactory, giving a maximum discrepancy of about 7%, for the case with the largest equilibrium flow, for which the asymptotic conditions (64) are more hardly satisfied.

For c β = 0.02, the agreement between numerical and analytical results is again very good but it emerges that, already in this β regime, the effect of the equilibrium flow is still stabilizing but quantitatively very modest, even for large flows. In particular, if compared to the results of Sec. III, Fig. 3 suggests that small ∆ ′ regimes are those where the influence of an out-of-plane flow on the linear stability of reconnecting modes can be most important.

V. CONCLUSIONS

We investigated the linear stability properties of inertial reconnection triggered by perturbing an equilibrium with a sheared magnetic field and an out-of-plane velocity jet.

Analytical dispersion relations have been derived in both the small and large ∆ ′ regimes.

Such dispersion relations, corresponding to Eq. ( 51) and (73), respectively, can be written, for finite β and neglecting terms of order d 2 e /d 2 i , as

π d e 1 -2 v a d 2 e d i γ + 1 + 2 v a d i 3/4 c β d 1/2 i (ak) 1/2 γ 1/2 -d β ak∆ ′ = 0, (75)
and

γ 3 = 2 π a 3 k 3 d e d 2 β G γ akd β c β (1 + 2vd i /a) 1/2 . ( 76
)
For the convenience of the reader, we recall that γ here indicates the growth rate, k the wave number, d i,e the ion and electron skin depth, c β = β/(1 + β), d β = d i c β and finally a and v, correspond to the amplitudes of the equilibrium magnetic and velocity fields, respectively.

The analysis of the above dispersion relations shows that increasing the jet amplitude results in reducing the growth rates, thus slowing the island evolution. The instability threshold, however, is not modified. The analytical predictions are confirmed, to a high degree of accuracy, by numerical simulations of the adopted four-field model.

The analysis lets us conclude that the stabilizing role of the out-of-plane flow is stronger in the small ∆ ′ regimes, where, for instance, sub-Alfvénic jets with a peak amplitude of about one tenth of the out-of-plane Alfvén speed are seen to cause a drop of the growth rate of about 20%. Also, if the jet has a characteristic thickness l smaller than the magnetic shear length, even slower equilibrium flow can affect the growth rate. In the large ∆ ′ regime, however, the influence of the velocity jet on the growth rate appears in general to be very weak, although still stabilizing. We remark also that, such stabilizing effect of out-of-plane differential flows is compatible with a similar behavior, described in Ref. [16] for a reduced MHD toroidal model.

On the other hand, if we restrict to the specific case of tokamak plasmas, the impact on the growth rate, of the class of equilibrium flows that we consider, appears to be negligible.

For instance, for a plasma with hydrogen ions, assuming B 0 = 30Kgauss, a = 0.1, v = 10Km/s and L = 5mm yields υd i = 0.014 ≪ 1. This indicates that collisionless tearing modes in tokamaks cannot be suppressed or effectively stabilized with the symmetric class of equilibrium flows that we considered, and suggests that asymmetric flows should be studied in this context. However, we do not exclude that the adopted symmetric equilibria might have a non-negligible impact in other situations in which out-of-plane flows and reconnection interact, such as for instance in the magnetotail [13,14], where Alfvénic or even super-Alfvénic flows are also observed.

The effect of the parallel compressibility on the growth rate for magnetostatic equilibria in the small ∆ ′ regime has also been analyzed. It emerges that, accounting for its presence in the dispersion relation, yields slightly lower growth rates, which fit better the numerical datas. In particular, this extends the dispersion relation of Ref. [18], which was independent on such effect.

We speculate that the inclusion of parallel compressibility also in the treatment of the large ∆ ′ dispersion relation, might help in reducing the analytical growth rates, which, as shown in Fig. 3, slightly exceed the numerical results.

We remark that, although applied to the Hamiltonian four-field model ( 1)-( 4), which describes inertial reconnection, our analysis can easily be extended to resistive tearing modes.

On the other hand, the inclusion of asymmetric equilibrium velocity profiles might require a substantial modification of the present linear analysis.

Finally, we anticipate that the natural question, about whether the presence of the out-ofplane velocity jet might influence also the nonlinear evolution of the reconnection process, will be the subject of a subsequent article and will complement the linear analysis described in the present paper.

  of the electrostatic potential φ and of the z (i.e. out-of-plane) component of the plasma velocity, denoted as v. All fields in the model have translational invariance along the z coordinate. The parameter d e corresponds to the electron skin depth, whereas the parameters c β and d β are defined as c β = β/(1 + β) and d β = d i c β , with d i indicating the ion skin depth. The canonical bracket is defined by [f, g] = ẑ • ∇f × ∇g for two generic functions f and g. The system (1)-(

FIG. 1 .

 1 FIG. 1. Comparison between values of growth rates obtained from the analytical dispersion relation (51) and from numerical simulations in the small ∆ ′ regime, for c β = 0.098 (top figure) and c β = 0.3 (bottom figure). The analytical results correspond to the solid line, whereas the numerical results are indicated by triangles. From the figure one can see that the numerical results are in quantitative agreement with the analytical predictions. In particular, the values of the growth rate decrease when increasing υd i , indicating the stabilizing role of the out-of-plane flow. The remaining values of the parameters for this figure are ∆ ′ = 0.9, d e = 0.2, d β = 0.39, a = 0.2. The corresponding numerical simulations have been carried out on the domain {(x, y) : -11.316 ≤ x ≤ 11.316, -1.25π ≤ y ≤ 1.25π}.

FIG. 3 .

 3 FIG. 3. Comparisons between values of growth rates obtained from the dispersion relation (73) and from numerical simulations in the large ∆ ′ regime. Squares (triangles) indicate the analytical (numerical) values. The top frame refers to the case c β = 0.1, whereas the bottom frame refers to c β = 0.02. The curves show a stabilizing effect due to the equilibrium flow, as in the small ∆ ′ regime. The effect, however, is quantitatively smaller and actually negligible for the case with the lowest value of c β . Parameters are ∆ ′ = 43.23, d e = 0.2, d β = 0.4, a = 0.2. The simulation domain is {(x, y) : -28.32 ≤ x ≤ 28.32, -30π ≤ y ≤ 30π}.

ACKNOWLEDGMENTS

The authors acknowledge useful discussions with Dario Borgogno, Cristel Chandre, Enzo Lazzaro, Fulvio Militello, Maurizio Ottaviani and François Waelbroeck. This work was supported by the European Community under the contracts of Association between EURATOM and ENEA and between EURATOM, CEA, and the French Research Federation for fusion studies. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Financial support was also received from the Agence Nationale de la Recherche (ANR GYPSI n. 2010 BLAN 941 03) and from the CNRS through the PEPS

IV. STABILITY ANALYSIS: LARGE ∆ ′ REGIME

For large ∆ ′ , the constant-ψ approximation adopted in the previous section is no longer valid. We resort then to different assumptions. In particular, our analysis mirrors that of Refs. [18,23]. In the latter references the reader can find a detailed explanation of the adopted procedure, so that here we describe only the main steps carried out to derive the dispersion relation.

The starting point is assuming the ordering

which allows to solve the linearized system ( 7)- (10)

A. Ideal MHD region

On scales over which d β and d e can be neglected, Eqs. ( 7) and ( 8) decouple from the system. Thus, as in the small ∆ ′ regime, one can find solutions for ψ and φ taking advantage of the assumption γ ≪ 1. For the large ∆ ′ case, it turns out to be convenient to consider in particular the solution for φ, which, in the vicinity of x = 0, using ( 22) and ( 24), reads

where φ 0 is an arbitrary constant and, again, we are considering the problem only on x ≥ 0.

B. Ion region

In the ion region the relevant system corresponds to that adopted for the outer region in the small ∆ ′ regime, and which we reformulate here in the following way:

which determines the solution for the out-of-plane perturbation of the magnetic field in the electron region.

In Fourier space, the matching with the solution in the ion region concerns the limit r → 0 of the electron region solution. In this limit the leading order solution of (70) [18] is given by

where Z 0 is an arbitrary constant and A e = Q -1 . On the other hand, the Fourier transform of the ion region solution (61) reads [18]

From ( 59) one also obtains that, for s → 0, V ∝ d 2 φ/dx 2 . The electron region solution (71), on the other hand, is also proportional to d 2 φ/dx 2 , because, from ( 68)-( 69), one obtains that Z ∝ p 2 φ. Consequently, the matching can be carried out directly between V from the ion region and Z from the electron region. In the limit 0 < ν ≪ 1, the result of the matching gives

Eq. ( 73) is the dispersion relation in the large ∆ ′ regime, including the effect of the equilibrium out-of-plane flow. For α = 1 one retrieves the result of Ref. [18], valid for static equilibria. As expected, in the limit of vanishing β, the effect of the out-of-plane equilibrium flow disappears and (73) reduces to the dispersion relation [30] g = 2 π

1/3

In order to show the effect of the flow on the growth rate and test the validity of the relation (73), a comparison between the analytical predictions of (73) and the growth rates obtained from numerical simulations is shown in Fig. The case with higher β (c β = 0.1) shows again a stabilizing effect of the equilibrium flow, as in the small ∆ ′ regime. We remark that, in the limit v → +∞, the dispersion relation (73)