The physics of a small-scale tearing mode in collisionless slab plasmas

Microtearing modes have been widely reported as a tearing parity electron temperature gradient driven plasma instability, which leads to fine scale tearing of the magnetic flux surfaces thereby resulting in reconnection of magnetic field lines and formation of magnetic islands. In slab geometry it has previously been shown that the drive mechanism requires a finite collision frequency. However, we find in linear gyrokinetic simulations that a collisionless fine-scale tearing parity instability exists even at low and zero collision frequency. Detailed studies reveal that these slab modes are also driven by electron temperature gradient but are sensitive to electron finite Larmor radius effects, and have a radial wavenumber much smaller than the binormal wavenumber, which is comparable to the ion Larmor radius. Furthermore, they exist even in the electrostatic limit and electromagnetic effects actually have a stabilising influence on this collisionless tearing mode. An analytic model shows that this collisionless small scale tearing mode is consistent with a tearing parity slab electron temperature gradient (ETG) mode, which can be more unstable than the twisting parity ETG mode that is often studied. This small-scale tearing parity mode can lead to magnetic islands, which, in turn, can influence turbulent transport in magnetised plasmas.


Introduction
In general, electromagnetic micro-instabilities in magnetised plasmas can be categorised as tearing or twisting parity modes. Tearing parity modes, in which the fluctuating parallel component of the magnetic potential is an even function about the rational surface, perturb the magnetic field to form magnetic islands. Twisting parity modes have a parallel component of the magnetic potential which is odd about the rational surface, and cause a rippling of the flux surface. Microtearing modes (MTMs) are a type of tearing parity micro-instability.
MTMs are a candidate for anomalous electron heat transport in tokamak plasmas [1,2,3]. They have been studied extensively since the 1960s. They are characterised by large toroidal and poloidal wavenumbers, comparable to the reciprocal of the ion Larmor radius. An early analytic linear model for MTMs was developed by J. F. Drake et. al [4]. In that work, the main drive mechanism is shown to come from the free energy in the electron temperature gradient, contributing to electromagnetic fluctuations. The collision frequency was shown to be a key factor in this drive mechanism. In the simple two dimensional sheared slab geometry, they studied the impact of collision frequency by dividing it into high-collisional, semi-collisional and collisionless regimes, and predicted MTMs to be stable at both low and high collision frequency. Numerical calculations by N. T. Gladd et. al [5] confirmed these slab results, demonstrating that a velocity dependent collision operator is essential for instability.
More recent studies [6,7] have observed MTMs in simulations neglecting the velocity dependence of the collision operator and even in the limit that the collision frequency tends to zero. Furthermore, gyrokinetic simulations have found microtearing modes can exist towards the edge of MAST tokamak plasmas and that these modes can be unstable at low collision frequency in toroidal geometry [7]. Gyrokinetic simulations [8,9,10] have also demonstrated unstable MTMs in the complete collisionless limit in a range of scenarios. These are at odds with the slab results presented in [5,4] but the mechanism is not yet fully understood. Understanding the collisionless drive mechanism is vital for clarifying the impact for transport in tokamak plasmas, especially those operating at higher temperature such as ITER.
In this article, we show that a fine scale collisionless tearing parity mode can, in fact, be unstable even in slab geometry. To identify the key physics, we develop an analytic model of this collisionless microtearing instability in slab geometry and conclude that the drive mechanism persists even in the electrostatic limit, with finite electron Larmor radius effects playing an important role. Probing the model in more detail, we show that the instability is the tearing parity branch of the electron temperature gradient (ETG) mode, which can be more unstable than the usual twisting parity ETG mode.
The paper is laid out as follows. In the next section we describe the magnetic geometry. In Section 3 we employ the GS2 code to demonstrate the existence of a microtearing instability in a collisionless slab. In Section 4 we discuss details of our analytic model, which we use to identify the main physics mechanisms in Section 5. We close in Section 6 with conclusions.

Slab geometry
It is convenient to define the slab geometry before our discussion of the physical plasma instability. We consider a simple infinite slab of plasma with magnetic field lines in the y − z plane and with density and temperature gradients in the x direction. The scale lengths are L −1 n = −dlnn/dx and L −1 T = −dlnT /dx, respectively. Using L T as the reference length, the normalised temperature gradient is defined as η = L n /L T . An external magnetic field B and a current density J are applied along the z direction, resulting in B = B 0 (ẑ + (x/L s )ŷ), where L s represents the scale length of the magnetic field shearing. We assume that |B y | |B z |, so restrict consideration to x L s .

GS2 simulation
GS2 is an initial value simulation code solving the gyrokinetic Vlasov-Maxwell equations using an implicit algorithm [11]. It employs local flux tubes and is designed to operate in a range of magnetic geometries including general tokamak, cylindrical and slab plasmas. We first employ GS2 (version v8.0.1 [12]) to benchmark the numerical results obtained by Gladd et al [5] in the slab geometry. Ion and electron temperatures are equal at the centre of the slab; however, the ion temperature gradient is zero while the electron temperature gradient is finite. There is also a finite density gradient and sheared magnetic field applied as described in Section 2. The scale lengths for the sheared magnetic field and for each species' temperature gradient and density gradient are L s , L T and L n , respectively. The mode frequency and growth rate in this paper are normalised to the electron diamagnetic frequency ω * e = k y v e ρ e /2L n . Here, the wavenumber k y and spatial coordinate x will be normalised to the ion gyro radius, As shown in figure (2) the linear GS2 results match well with Gladd's model [5] in the collisional regime, and both demonstrate the drive from the electron temperature gradient. However, in the very low collision frequency regime, GS2 reveals an unexpected tearing parity instability in this slab geometry. The frequency is not continuous between the collisional and collisionless regions, indicating that they are different instability branches. The real and imaginary parts of the normalised electrostatic potentialφ and normalised parallel magnetic potentialĀ are shown in figure (3).Ā has an even symmetry whileφ is odd, which is a defining feature of tearing parity modes. The collisionless one has a narrower structure, thus the characteristic radial wavenumber k x ρ i is much larger for the collisionless branch than for the collisional one (see figure (4)). In fact, in these GS2 simulations we note that capturing the collisionless instability requires challenging numerical settings. The parallel grid extent and resolution needs to be sufficiently high to capture the unstable mode accurately. In GS2, the parallel flux tube extent is controlled by nperiod, while ntheta defines the grid resolution within each 2π period. In our simulations, the collisionless branch requires nperiod = 128 and ntheta = 8, while the collisional branch is well converged for nperiod = 32 and ntheta = 8.
Note that the condition k x ρ i 1 is assumed in the derivation of references [4,5]. This enables a Gamma function expansion in the quasi-neutrality equation, ignoring the finite Larmor radius effects from electrons. Figure (4) tests the validity of this assumption for the range of collision frequencies, with η e = 5.0 and η e = 7.0. It shows that electron temperature gradient η e has very little influence on the value of k x ρ i but collision frequency ν/ω * e has a big impact. In the collisional regime, k x ρ i remains small and the approximation is valid; however, this is not the case for the collisionless regime. Specifically, in the collisionless regime, k x ρ i is about 40 times larger, which leads to k x ρ e approaching 1. This gives the first clue that the finite Larmor radius effects from electrons might be an important factor for the collisionless mode seen here. This finding informs a new reduced gyrokinetic model, to be derived in Section 4.
To further test this point, we examined the influence of electron finite Larmor radius effects directly by probing a Bessel function parameter in GS2. This parameter, α, enables a suppression of finite Larmor radius effects in the gyro-averaging Bessel function J 0 (αk). α = 1 is the default case capturing full gyrokinetic physics. Turning α down towards zero is equivalent to turning off the finite Larmor radius effects for the given species in GS2 simulations. Figure (5) shows the effects of the Bessel factor on the collisionless and collisional branches in GS2 simulation. Here the electron temperature gradient is η e = 5.0 and the collision frequency is set to ν/ω * e = 0.1 for the collisionless simulation. Note that the gyro-averaging provides velocity dependent dissipation, which is in some sense similar to the collision operator. Whilst the collisionality is an important factor in the collisional slab model, this indicates that for the collisionless instability, the electron finite Larmor radius effects are required to confine the mode; these are neglected in the collisional model. On the contrary, the collisional instability is insensitive to the electron finite Larmor radius effects. Focusing on the zero collision frequency limit and varying the ratio of plasma pressure to magnetic pressure, β, we found that this mode persists in the electrostatic limit, as shown in figure (6). Indeed, it is more unstable at lower β, and still exists when  β = 0. This means that electromagnetic effects are stabilising for this mode, which is fundamentally electrostatic in nature. Meanwhile, GS2 simulations with kinetic and adiabatic ions demonstrate that the ion treatment has little impact on the collisionless mode. This provides conclusive evidence that the main drive comes from electrostatic

Modelling in slab geometry
In order to provide a physics interpretation of the GS2 simulation at low collision frequency, we derive eigenmode equations valid in this limit. Here we present two models derived from gyrokinetic theory. Section 4.1 describes a simple case, focusing on zero collision frequency and zero β, which demonstrates just the essential physics. Section 4.2 considers more complicated factors including finite but small collision frequency and electromagnetic effects, which can be compared in more detail with GS2 results and help give a good foundation for future work.

Electrostatic model at zero collision frequency
Informed by the earlier GS2 results, we adopt an adiabatic ion response but treat electrons kinetically, retaining finite Larmor radius effects. We first consider the electrostatic limit, with perturbations only in the electrostatic potential φ. In Fourier space the gyrokinetic equation for electrons yields: Assuming small k ⊥ ρ ⊥ , expanding Bessel functions to second order and conducting an inverse Fourier transform to real space, equations (1) and (2) become: Substituting (3) into (4) and normalising the variables asω = ω/ω * e ,k y = k y ρ e , x = x/ρ e ,φ = eφ/T , we have in which s = v /v e , t = v ⊥ /v e , = L n /L s and η = L n /L T . Please note that we have normalised lengths to the electron Larmor radius rather than the ion Larmor radius in the previous sections. We considerk y 1 in which case it can be neglected. Neglecting third and fourth orders of the expansion in k x ρ e , simplification of this equation yields a second order differential equation for the electrostatic potential of the form C 0φ + C 1φ + C 2φ = 0, where primes denote the differential with respect tox and the coefficients C 0 , C 1 and C 2 are Here, Z m,n = Z m,n (ω/2 x) is a generalised plasma dispersion function It can be shown that Z 0,0 (α) = iπW(α) where W(α) is the Faddeeva function. When mn = 0, there is a pair of recurrence relations which can be used to relate Z m,n (α) to Z 0,0 (α): When near the centre of the slab, wherex = 0, the above coefficients are welldefined with the limit of The forms of these coefficients show that, when normalising to the diamagnetic frequency ω * e , the mode frequency and growth rate are mostly sensitive to magnetic shear scale length and electron temperature gradient. Numerical solutions of this second order differential equation forφ are presented in Section 5.

Electromagnetic model with finite Lorentz collision operator
The above electrostatic model is valid only when the collision frequency is zero. To compare with the GS2 results along the low collision frequency range, we consider a classic Lorentz collision operator consisting of pitch-angle scattering C(ν) = − iν where ν is the collision frequency and ξ = v /v is the pitch angle. Furthermore, it is useful to explore the influence of electromagnetic effects including β, which will show the tendency of this mode to form magnetic islands. Before we address this, we first update the above model. When including the parallel magnetic potential and the Lorentz collision operator, the gyrokinetic equation can be rewritten as Again, expanding the Bessel function and conducting an inverse Fourier transform results in Note that the perturbation of the electron distribution function depends on both space and velocity g(x) = g(x, v, ξ). Expanding the distribution function in an orthogonal polynomial series, g(x, v, ξ) = ∞ n=0 h n (x, v)P n (ξ) in which P n is the Legendre polynomial of order n, we have Applying the orthogonality relations for Legendre polynomials and integrating over pitch angle ξ from −1 to 1 on both sides yields a set of equations To derive a tractable model from the above, we adopt a matrix approach. The equation (18) can be written in the matrix form as M · h = D. The two dimensional matrix M is an infinite tridiagonal matrix, of which the n-th row (n starts from 1) is The column vector h starts from the h 0 term and the column vector D represents the driving terms in the right hand side of equation (18). Note that the main difference between this model and collisional model in reference [5] lays in D. Without finite Larmor radius effects, D becomes a scalar thus it is possible to present h terms in a continued fraction as in their model. In our model, re-writing as and noting that d 2 M /dx 2 = 0, we have Substituting the Legendre series of g(x, v, ξ) into the quasi-neutrality equation (2) and expanding the Bessel function as before, we have We define two row vectors as Therefore the quasi-neutrality equation (21) becomes To complete the set of electromagnetic equations, Ampere's Law is also needed. The parallel perturbed current density in Fourier space iŝ Similarly, we expand the Bessel function, define another two row vectors as and follow the same procedure above. The parallel perturbed current density in real space becomes Therefore, Substituting the expressions for h in equation (19) and d 2 h/dx 2 in equation (20) into equations (24) and (29), yields the final expressions incorporating both electromagnetic effects and a Lorentz collision operator.
Simplifying the parameters in the drive term D and normalising the results using the same parameters as in the electrostatic model above, as well as u = v/v e and A = A /ρ e B. The normalised matrix equations becomē These represent an infinite tridiagonal matrix and a column vector, respectively. In practice, they will be truncated, albeit at a large size during the calculation. Equations (30) and (31) will eventually lead to a system of two simultaneous second order differential equations for both electrostatic potential φ and parallel magnetic potential A . Generally, however, the coefficients for each of the terms are not easy to simplify and reveal the insight of physics, except in some special cases. For the collisional microtearing theory where the finite Larmor radius effects are not considered, the drive term D has non-zero elements only in the first two terms; thus only the left top four elements in the inverse matrix of M will contribute to the results. In this case, the matrix products can be simplified to produce the continued fraction and electron parallel conductivity in reference [5]. In another simple case, when collision frequency ν = 0, at the centre plane of the slab x = 0, the tridiagonal matrix M becomes diagonal, and the calculation will be significantly simplified. Dropping electromagnetic terms, the equation becomes equivalent to the electrostatic model discussed in the previous section.

Numerical results and discussion
We have established two reduced models for the collisionless micro-scale tearing instability considering finite Larmor radius effects from electrons. Based on the electrostatic eigenmode equations (6), (7) and (8), and electromagnetic eigenmode equations (30) and (31), we have developed two codes to calculate the complex eigenmode frequency ω for each of these models. The algorithm of the codes is an eigensolver based on an iteration method described in Chapter 5.4 of reference [13]. Here, for our electromagnetic equations with our chosen parameters, the complex frequency ω is typically converged to a relative tolerance of 10 −6 when the matrices are of size 30 by 30. The boundary conditions for the tearing instability for both models are spatially localised tearing parity requirements |φ|, |A | → 0 as |x| → ∞ , φ(0) = 0 , dA dx x=0 = 0 (40) Please note that we did not assume any parity in the derivation; thus our model is also capable of looking for twisting parity solutions, whose boundary conditions are For the collisionless limit ν = 0, solving either electrostatic or electromagnetic eigenmode equations gives results close to those of GS2 (see figure (8) for the point at ν = 0). We know that the electrostatic model describes the electron temperature gradient (ETG) mode. However, the ETG mode is usually considered to be a twisting parity mode. Nevertheless, as with any eigenmode problem, there is a family of solutions (harmonics) with alternating parity, in which the twisting parity is the fundamental harmonic. In fact, eigenmodes of different harmonics can co-exist and there is no physical reason why the fundamental one should be the most unstable. Indeed, such phenomena, where the higher harmonics are more unstable, were found previously in both ion temperature gradient (ITG) mode [14,15] and ETG mode calculations [16,17]; studies [18,19] reported the existence of unstable high order ballooning (twisting parity) modes and that parity transition can happen under certain parameters; studies [20,21] demonstrated the excitation of higher harmonics and the parity mixture under certain scenarios. For our electrostatic model, figure (7) shows different eigenvalues of which the even and odd harmonics are twisting and tearing modes, respectively, solved with corresponding boundary conditions. Our results show that the most unstable mode in this case is the third order harmonic, which corresponds to the collisionless tearing parity instability we have found. When including the collision frequency, the results from our electromagnetic model are similar to those obtained from GS2, as shown in figure (8). Both models show that the this mode is driven by electron temperature gradient, consistent with the identification as an ETG. The mode growth rate decreases as the collision frequency rises, but the growth rate in GS2 has a stronger variation and switches to a different harmonic at ν = 1.0 for η e = 5.0 and ν = 0.7 for η e = 7.0, as indicated by the jump in frequency. Note that our electromagnetic model results are consistent with the third order harmonic in figure (7). One reason for the difference between our model and GS2 might be that the collision operator in GS2 differs from our model [22,23]. Though quantitatively slightly different, both our model and GS2 show that this collisionless tearing parity instability tends to be stabilised by collisions. To summarise, we conclude that the collisionless tearing parity instability found here in slab geometry is a form of ETG instability, with a different drive mechanism to the standard collisional slab MTM.  Although the underlying mechanism behind the collisionless tearing parity instability studied here is different from the collisional MTM, it still leads to magnetic reconnection and the formation of magnetic islands. Both GS2 and our electromagnetic model provide the mode structure for the collisionless instability as shown in figure (3). We can calculate the flux surfaces of the magnetic field from the magnetic potential A . The structure of magnetic field lines is given by contours of the flux The contour plot of constant levels of ψ(x, y) gives the flux surfaces and hence the magnetic structure. Figure (9) shows the island structure for the collisionless mode when ν = 0. However, note that the amplitude of A is arbitrary in our linear model, so the width of the island is not determined. From left to right, the three panels of figure (9) are examples to show that under the same parameters the island shape can become more contorted as its size grows from the order of electron Larmor radius ρ e to ion Larmor radius ρ i . It can also be found that as the island width grows, a secondary island arises in the vicinity of the X-point. We believe that the inflection points of A will finally provide a limit for the maximum island width. How this will affect the particle and heat transport is to be answered in future work.

Conclusion
We have shown that there is a collisionless micro-scale tearing parity instability that can drive reconnection even in the absence of collisions. We have established two models considering electron finite Larmor radius effects to interpret the physics of this mode, which is shown to be stabilised when the collision frequency increases. We identify the collisionless mode as a tearing parity harmonic of the conventional slab ETG mode, which is the most unstable harmonic for our parameters. The electromagnetic component results in magnetic islands. Our result stands as an example to show that tearing parity modes can arise from a whole range of different drives, and there may be other possible ways to get small scale tearing parity modes. These can have an impact on the transport and can be very challenging to resolve numerically, posing problems for attempts to simulate them. On the other hand, even if the tearing parity eigenmodes are not the most unstable harmonic linearly, it still may be possible that tearing harmonics can play a role nonlinearly, leading to a background degradation to the confining magnetic field everywhere that such instabilities exist.
The remaining questions in our research include why and in what parameter range does the tearing harmonic become the most unstable ETG mode in the electrostatic model, and how does this collisionless mode behave in toroidal geometry. The electromagnetic model (30) and (31) we obtain is a rather complicated expression. Further simplification may reveal more physical insight. This work is funded by China Scholarship Council and the University of York. It is also part of TDoTP project funded by EPSRC (EP/R034737/1).