The linear tearing instability in three dimensional, toroidal gyrokinetic simulations
William Hornsby, Pierluigi Migliano, Rico Bucholz, Lukas Kroenert, Arthur Peeters, David Zarzoso, Emanuele Poli, Francis Casson
TThe linear tearing instability in three dimensional, toroidal gyrokinetic simulations.
W.A. Hornsby, P. Migliano, R. Bucholz, L. Kroenert, A.G. Peeters
Theoretical Physics V, Dept. of Physics, Universitaet Bayreuth, Bayreuth, Germany, D-95447 ∗ D.Zarzoso, E. Poli
Max-Planck-Institut f¨ur Plasmaphysik, Boltzmannstrasse 2, D-85748 Garching bei M¨unchen, Germany
F.J. Casson
CCFE, Culham Science Centre, Abingdon,Oxon, OX14 3DB, UK (Dated: October 30, 2018)Linear gyro-kinetic simulations of the classical tearing mode in three-dimensional toroidal ge-ometry were performed using the global gyro kinetic turbulence code, GKW . The results werebenchmarked against a cylindrical ideal MHD and analytical theory calculations. The stability,growth rate and frequency of the mode were investigated by varying the current profile, collision-ality and the pressure gradients. Both collision-less and semi-collisional tearing modes were foundwith a smooth transition between the two. A residual, finite, rotation frequency of the mode evenin the absense of a pressure gradient is observed which is attributed to toroidal finite Larmor-radiuseffects. When a pressure gradient is present at low collisionality, the mode rotates at the expectedelectron diamagnetic frequency. However the island rotation reverses direction at high collisionality.The growth rate is found to follow a η / scaling with collisional resistivity in the semi-collisionalregime, closely following the semi-collisional scaling found by Fitzpatrick. The stability of the modeclosely follows the stability using resistive MHD theory, however a modification due to toroidalcoupling and pressure effects is seen. ∗ [email protected] a r X i v : . [ phy s i c s . p l a s m - ph ] A ug FIG. 1. Contour plots of (left) the parallel-radial mode structure of A (cid:107) and (right) the electrostatic potential, φ for a m = 2, n = 1tearing mode. The vertical dashed line denotes the position of the q = 2 rational surface. I. INTRODUCTION
A current density profile in a plasma can lead to instabilities that are able to change the topology of the magneticfield via the process of magnetic reconnection [1]. Resistive instabilities, and in particular tearing instabilities [2–4]are found in a wide variety of astrophysical and laboratory situations [5].In a magnetically confined tokamak plasma, reconnecting instabilities can form magnetic islands. This has adetrimental effect on the plasma confinement caused by the introduction of a radial component of the magnetic field[6]. The tearing mode can trigger a disruption, something which would be catastrophic to the ITER tokamak andtherefore the tearing mode, and specifically, the neoclassical tearing mode (NTM) [7, 8] is expected to be the limitingfactor on its operational plasma β [10].Much work has been performed numerically in studying tearing mode stability in both cylindrical and toroidalgeometry [11], using the Magnetohydrodynamic (MHD) framework. However for large current and future tokamaksa fully kinetic description of the physics in the resonant layer is needed for both the linear and non-linear regimes.Linear theory is generally considered to be valid when the generated magnetic island has a width that is smallerthan the singular layer. Above this width, nonlinear effects slow the growth of the island to an algebraic scaling [12]and eventually cause the island to saturate [13]. Recent work studying the effects of electromagnetic turbulence onthe evolution of the tearing mode indicates a turbulent modification of the resonant layer, extending the suitability oflinear theory to larger island sizes[14], while it has also been seen that a magnetic island can have a profound effecton turbulence and transport [15, 16].Within the singular layer (Parallel wave-vector of the mode, k || = 0) the assumptions of ideal MHD break downand magnetic reconnection may take place. For a collisionless mode the mechanism of reconnection is the electroninertia, where the singular layer width, and in turn the growth rate, are related closely to the electron skin depth [17].When collisions become significant, it is plasma resistivity which determines the singular layer width, producing aclassical collisional tearing [3] or, more relevant for present day tokamaks, a semi-collisional mode [18, 19]. The modefrequency, and in particular, its direction is vital in predicting the nonlinear stabity of the NTM. Specifically the signof the mode frequency, whether the mode rotates in the ion or electron diamagnetic direction, determines whetherthe polarisation current has a stabilising or destabilising effect [9, 23] particularly when the magnetic island is small[24].In this work the linear mode stability, growth rate and frequency is studied for the first time in three-dimensionaltoroidal geometry using the gyro-kinetic framework of equations. Parameters appropriate to modern and futureTokamak experiments are used. Gyrokinetics has been highly successful when applied to the study of drift waves,turbulence and transport but its use in the study of large scale MHD instabilities is still in its infancy. Recent work hasstudied two dimensional magnetic reconnection in a highly magnetised system [25–27] and the ideal kink instability[28].The paper is organised as follows. In Sec.II we describe the equations solved, while in Sec. III we outline howthe tearing mode is implemented within the gyrokinetic framework. Sec. IV outlines the parameters used, and theresolution requirements. Sec. V and VI presents the scaling of the tearing mode growth rate and mode frequency withcollisionality, and with the variation of the background density, temperature and current profiles. II. GYROKINETIC FRAMEWORK
The self-consistent treatment of the tearing mode requires a radial profile in geometry and thermodynamic quan-tities. We use here the global version of the gyrokinetic code, GKW [29]. This is advancement on the already widelyused flux-tube variant [30].The code solves the gyrokinetic set of equations. The full details can be found in [30] and references found therein.Here we outline the basic set of equations that are solved and in the next section the modification and assumptionsused in driving the tearing instability.The delta- f approximation is used. The distribution function is split into a background F and a perturbed dis-tribution f . The final equation for the perturbed distribution function f , for each species, s , can be written in theform ∂g∂t + ( v (cid:107) b + v D ) · ∇ f − µBm B · ∇ BB ∂f∂v (cid:107) = S, (1)where S is the source term which is determined by the analytically known, background distribution function, µ isthe magnetic moment, v || is the velocity along the magnetic field, B is the magnetic field strength, m and Z are theparticle mass and charge number respectively. Here, g = f + ( Ze/T ) v (cid:107) (cid:104) A (cid:107) (cid:105) F M is used to absorb the time derivativeof the parallel vector potential ∂A (cid:107) /∂t which enters the equations through Amp`eres law.The velocities in Eq. (1) are from left to right: the parallel motion along the unperturbed field ( v (cid:107) b ), the drift motiondue to the inhomogeneous field ( v D = Ze (cid:20) mv (cid:107) B + µ (cid:21) B ×∇ BB ), and the motion due to the perturbed electromagneticfield ( v χ = b ×∇ χB , where χ = (cid:104) φ (cid:105) − v || (cid:104) A || (cid:105) ). The latter is the combination of the E × B velocity ( v E = b × ∇(cid:104) φ (cid:105) /B )and the parallel motion along the perturbed field line ( v δB = − b × ∇ v (cid:107) (cid:104) A (cid:107) (cid:105) /B ). Here, the angled brackets denotegyro-averaged quantities.The background is assumed to be a shifted Maxwellian ( F M ), with particle density ( n ( r )) and temperature ( T ( r ))and toroidal flow velocity ( ω φ ( r )) F Ms = n s π / v exp (cid:20) − ( v (cid:107) − ( RB t /B ) ω φ ) + 2 µB/mv (cid:21) , (2)which determines the source term, S = − ( v χ ) · (cid:20) ∇ nn + (cid:18) v (cid:107) v + ( µB ) T − (cid:19) ∇ TT (3)+ mv (cid:107) T RB t B ∇ ω φ (cid:21) F M − ZeT [ v (cid:107) b + v D ] · ∇(cid:104) φ (cid:105) F M . (4)The thermal velocity v th ≡ (cid:112) T /m , and the major radius ( R ) are use to normalise the length and time scales. Usingstandard gyro-kinetic ordering, the length scale of perturbations along the field line ( R ∇ (cid:107) ≈
1) are significantly longerthan those perpendicular to the field ( R ∇ ⊥ ≈ /ρ ∗ ). Here, ρ ∗ = ρ i /R is the normalised ion Larmor radius (where ρ i = m i v th /eB and v th = (cid:112) T i /m i ).The electrostatic potential is calculated from the gyro-kinetic quasineutrality equation, (cid:88) s Z s e (cid:90) (cid:104) f (cid:105) † d v + (cid:88) s Z s e T s (cid:90) ( (cid:104)(cid:104) φ (cid:105)(cid:105) † − φ ) F Ms ( v ) d v = 0 (5)where the first term represents the perturbed charge density and the second the polarisation density, which is onlycalculated from the local Maxwellian. The dagger operator represents the conjugate gyro-average operator. Similarlythe parallel vector potential is calculated via Amp`eres law, − ∇ A (cid:107) + (cid:88) s µ Z s e T s (cid:90) d v v (cid:107) (cid:104)(cid:104) A (cid:107) (cid:105)(cid:105) † F Ms = (cid:88) s µ Z s eT s (cid:90) d v v (cid:107) (cid:104) g s (cid:105) † (6)The gyro-average is then calculated as a numerical average over a ring with a fixed radius equal to the Larmorradius using a 32 point interpolation defined by, (cid:104) G (cid:105) ( X ) = 12 π (cid:73) d α G ( X + ρ ) (7)where α is the gyro-angle, X is the position of the gyro-center and ρ is the gyro-radius. This gyro-average is used inboth the evolution equation of the distribution function, as well as in the quasineutrality and Amp´ere equations andon both the electron and ion species. The nonlinearity ( v χ · ∇ g ) is neglected in this study and its effects on magneticisland evolution and saturation will be left for a future publication.GKW uses straight field line Hamada [32] coordinates ( s, ζ, ψ ) where s is the coordinate along the magnetic fieldand ζ is the generalised toroidal angle and ψ = r/R ref is the normalised radial coordinate. For circular concentricsurfaces [31], the transformation of poloidal and toroidal angle to these coordinates is given by [30] ( s, ζ ) = ( θ/ π, [ qθ − φ ] / π ). The wave vector of the island is k Iζ ρ i = 2 πnρ ∗ , where n is the toroidal mode number. GKW uses a Fourierrepresentation in the binormal direction, perpendicular to the magnetic field. The radial and parallel directions aretreated using a high order finite difference scheme so as to include global profile effects utilising Dirichelet boundaryconditions in the radial direction, and open field line boundaries in the parallel direction. An example of the radialprofiles used is seen in the middle panel of Fig. 2.The collision operator is particularly important to the growth and development of the tearing mode, as it can beresistive in nature. In GKW the full linearised Landau collision operator with momentum and energy conservingterms is implemented. However here, unless stated otherwise, we utilise only the pitch-angle scattering part. Thedifferential part of the operator has the following form [33] C ( f a ) = (cid:88) b v sin θ ∂∂θ (cid:20) sin θD a/bθθ v ∂f a ∂θ (cid:21) . (8)Here the sum is over all species b and θ denotes the particle pitch-angle ( v || = v cos θ ).The coefficients can be obtained from the literature [33] and written in normalised form, suppressing the subscriptN for the normalised quantities: D a/bθθ = (cid:88) b Γ a/b v a (cid:20)(cid:18) − v b (cid:19) erf( v b ) + 1 v b erf (cid:48) ( v b ) (cid:21) , (9)where erf and erf (cid:48) are the standard definition of the error function and its derivative. Finally, the normalised collisionfrequency, Γ a/b is given byΓ a/b = R ref n b Z a Z b e ln Λ a/b π(cid:15) m a v tha , n b is the scattering species number density, Z is the relative charge of the species andln Λ is the Coulomb logarithm. v a and v b are the velocities of the scattered and scattering species respectively. III. TEARING MODE IMPLEMENTATION
The tearing mode instability is driven by a non-homogeneous current density. In our implementation this isintroduced by applying an electron flow profile in the equilibrium. The electron flow velocity is calculated self-consistently from the imposed q-profile analogous to the method used in [18] for a kinetic calculation in slab geometry.We assume that the background current is carried by the electrons, J = − n e v e e , i.e. that the electron flow velocity v e is larger than the ion flow and that it is a flux function. The current gradient is therefore related to the gradientin the electron flow, ∂u e /∂ψ , this profile is introduced to the code where ∇ ω φ is written in equation 4.We assume v e to be significantly smaller than the electron thermal velocity, thus the only term in Eq. 4 of interestis, ∇ F Me = − v || v th ∇ u e F Me . | A || | q , s , R / L n −3 ∂ A || / ∂ ψ GKW psil=0.03GKW psil=0.0MHDGKW, ε =1/6 ψ [r/R] | φ | / n i Shear, (1/q)(dq/d ψ )Safety factor, qR/L n | φ ||n i | FIG. 2. The radial profiles of (Top) the parallel vector potential as calculated by GKW and the ideal MHD solution incylindrical geometry calculated using the shooting method. The shift in the eigenmode is compatible with previous resultswhen comparing the cyclindrical to toroidal eigenfunction.[11] (Middle top) The safety factor (q),magnetic shear (ˆ s ) and, whenpresent, density gradient ( R/L n ) profiles (Middle bottom) the radial second derivative, ∂ A || /∂ψ , indicating the position ofthe resistive layer. (Bottom) are the electrostatic potential ( φ ) and perturbed ion density ( n i ) profiles. The electron profile is related to the q-profile via Amp´eres law, (cid:82) (cid:126)B · ds = µ (cid:82) jd A , which in circular geometryand assuming a low aspect ratio, becomes,2 πrµ B p F ( r ) = 2 πµ (cid:90) r r (cid:48) dr (cid:48) j ( r (cid:48) ) , (10)where the function F ( r ) = 1 / π (cid:82) dθ/ (1 + (cid:15) cos θ ) is approximated to be close to unity and subscript zeros representconstant values. We have used the definition q = π (cid:82) rB t RB p = rB t R B p F ( r ) and therefore the poloidal magnetic field, B p = B t R F ( r ) rq . The gradient of the background current is given by the expression, djdψ = B t µ R ddψ (cid:20) ψ ddψ (cid:20) ψ q (cid:21)(cid:21) (11)In this paper we make use of the Wesson analytic form of the current profile [36] which is defined by the expression, j = j (cid:32) − ψ (cid:18) Ra (cid:19) (cid:33) ν (12)which, in turn, gives a q -profile of the form, q = q a r /a − (1 − r /a ) ν +1 (13)where q a is the safety factor at ψ = a/R . This is related to the q on the axis, q , by q a /q = ν + 1, where ν is aninteger that determines that peaking of the current gradient.The density and temperature profiles have the radial form, ∂n s ∂ψ = 12 RL ns (tanh ( x − x + ∆ x ) /w − tanh ( x − x + ∆ x ) /w ) ∂T s ∂ψ = 12 RL T s (tanh ( x − x + ∆ x ) /w − tanh ( x − x + ∆ x ) /w )where R/L n = − ( R/n ) ∂n/∂ψ and R/L T = − ( R/T ) ∂T /∂ψ are the logarithmic density and temperature gradients atthe reference radius, x , ∆ x is the profile width and w is the rising width. IV. PARAMETERS
In this section we study the linear mode structure, a two-dimensional example of which is plotted in Fig. 1, showingthe radially elognated eigenfunction of the parallel vector potential (left panel) and the electrostatic potential ( φ )which is highly localised around the singular layer at the q = 2 rational surface (dashed line). The growth rates andmode rotation frequencies are calculated for the following parameters unless stated otherwise in the text:Simulation parameters.Aspect ratio, R/a 3Electron beta, β e − (0 . ν m i /m e q a k ζ ρ i ρ ∗ = ρ i /R T i = T e max( v || /µ ) 4/8Care must be taken in resolving the relevant scale lengths in each of the domain directions. The number of gridpoints in the poloidal, parallel velocity and µ directions were, N s = 32, N vp = 64, N µ = 16 respectively which werefound to give converged growth rates to within a couple of percent. A single toroidal mode is considered, namely withthe n = 1 mode number, the number of resonant poloidal modes within the domain is determined by the q -profile( q = m/n ); as such the possibilty of mode coupling effects remain, particularly the toroidal coupling of modes withvarying poloidal mode numbers (m) which is known to have a destabilizing effect [38]. The radial resolution requiredis dependent on the physics of the resonant layer. In the collisionless tearing mode the relevant scale length is theelectron skin depth, δ e = c/ω pe , where ω pe = (cid:112) n e e /m e (cid:15) is the electron plasma frequency and c is the speed of light.The skin depth is, δ e a = ρ ∗ (cid:113) mime √ β e Ra . The β e is the electron β defined by β e = n e T e / ( B / µ ), where n e and T e areequilibrium density and temperature and B is the magnetic field strength at the outboard midplane. For collisionalmodes, the resistivity defines the width of the layer. Unless otherwise stated, we have used, N x = 512 radial gridpoints between ψ l = 0 .
03 and ψ h = 0 . ρ ∗ = 0 . ∼ / β of 0 . (cid:15) → R/a = 10, closer to the cylindrical limit, whichshows that the eigenfunction tends toward this cylindrical form as expected. The radial density and geometry profilesused are shown in middle panel. The third panel of Fig. 2 shows the second derivative of the radial profile of thevector potential. This shows clearly the position of the singular layer, the position where ideal MHD breaks down a q a / q stableq=1q=2 no resonant surface q =1 A || , φ [ a r b ] S a f e t y f a c t o r , q −3 ψ [ ε /R] ∂ A || / ∂ ψ A || φ ∂ A || / ∂ x ∂ A || / ∂ x FIG. 3. (top) Stability of the classical tearing mode varying the current peaking parameter, ν = q a /q − q a , normalised collisionality, ν ∗ = 1 .
0. The lines correspond to the m = 2 , n = 1 stability boundaries asfound by Wesson et.al [36]. Solid circles denote an unstable mode, while crosses denote stability. (bottom) The radial modeprofiles for a coupled mode with q a = 3 . ν = 1. Inlaid in the top panel are the profiles of the first (dashed line) andsecond (solid line) radial derivatives of A || close to the q = 2 rational surface. and a discontinuity in the derivative is produced. The bottom panel also shows the radial profile of the electrostaticpotential ( φ ), showing the highly localised electric field at the rational surface, and the perturbed ion density profile( n i ). V. TEARING MODE STABILITY.
The upper panel of Fig. 3 shows the stability of the tearing mode as a function of the current peaking parameter, ν and the safety factor at the plasma edge, q a . This is analogous to the diagram shown in [36, 42], which shows thestability of the n = 1, m = 2 mode in cyclindrical geometry, n denoting the toroidal mode number, while m is thepoloidal mode number.The analysis in [36] considered a single m = 2, n = 1 mode, the outlines of which are also plotted in overlay (blacklines). However we are unable to decompose the modes in this manner and consider, for example a m = 3, n = 1mode in isolation. The different coloured points denote the dominant mode within a particular simulation determinedby the amplitude of the mode at the relevant rational surface.It is evident that there is good agreement with previous work, particularly when there exists only an m = 2, n = 1mode within the domain. However, the results are complicated by the toroidal coupling of modes. Toroidal couplinghas the effect of destabilising modes that would otherwise be stable [38, 39], for example in lower panel of Fig. 3 we seethe radial eigenstructure of a combined m = 2 , n = 1 and a m = 3 , n = 1 mode, the latter of which would otherwisebe stable. The point at q a = 4, ν = 2 has only a q = 3 rational surface within the domain and is found to be stable.This is further manifested in the growth rate where this is plotted if Fig. 4 for the combined m = 2 , m = 2 , n = 1 mode.The growth rate of the tearing mode is a function of the well known parameter [3],∆ (cid:48) = 1 A || ∂A || ∂r (cid:12)(cid:12)(cid:12)(cid:12) r + s r − s (14)across the singular layer, whose position is denoted by r s and upper and lower boundaries by r + s and r − s . Thisparameter represents the discontinuity in ∂A || /∂r and measures whether mode growth is energetically favourable(∆ (cid:48) > (cid:48) , measuring the growthrate allows us to study the mode stability. As such the above difference in growth rate between the coupled tearingmode and the pure m = 2 , n = 1 shows that the toroidal coupling in this case has a destabilising effect (Increasing∆ (cid:48) ), increasing the growth rate by a factor of 2.It can be seen that points along the line q = 1 are stable to the tearing mode. For the collisionality used inthe figure ( ν ∗ ∼ . S = 5 . ), a toroidal MHD analysis (Fig.8 in [42]) predicts that points along this line wouldindeed be stable, and the results from GKW are consistent with this, however, the mode is stabilised at a much lowerLundqvist number, S than is predicted in [42]. S is defined as the ratio of the Alven time ( τ A = a √ µ n i m i /B ) andthe resistive time scale τ R = µ a /η , where η = m e ν ei / ( n e e ). VI. COLLISIONALITY EFFECTS ON GROWTH RATE AND ROTATION.
Fig. 4 shows the growth rate of the tearing mode as a function of the normalised collisionality (normalised to thetrapping-detrapping rate, ν N = ν ei q/(cid:15) / ), corresponding to a Ludquist number, S, range between 10 and 2 × .A scan with no background density gradient (Black circles) and with a background gradient scale length of R/L n =1 . n x = 1024), the growth rate changed by less than 3%from the value using n x = 512.As the collisionality is increased (and as such the resistivity) we see the expected increase in the mode growth rate.The (blue dashed) line represents the condition γ = ν ei and marks the boundary between the collisionless and semi-collisional mode regimes, at this point collisions limit the electron response, which allows them to be accelerated, andbroadens the resonant layer width. The semi-collisional theory only holds in the limit when γ < ν ei and the transitionbetween semi-collisional and collisionless modes is clearly seen in our simulations. The growth rate obtained deviatessignificantly from the collisionless rate once the growth rate is greater than the electron-ion collision frequency.Plotted are the analytic scalings for the resistive (Green dashes, η / ) and the semi-collisional tearing mode (Blackdashes, η / ) obtained using a kinetic treatment by Drake and Lee [18]. However both of these give significantlystronger scalings than our calculated values. Good agreement, particularly when a density gradient is present, isfound when we consider the semi-collisional scaling at high ∆ (cid:48) as given by the three fluid, two dimensional slabcalculation by Fitzpatrick [19]. This predicts a scaling of the growth rate with η / (Light grey dashes). Thiscorresponds to regime III in Fig.1 of [19] and was also found in a semi-collisional study of m = 1 modes [20]. Thevisco-resistive scaling as calculated by Porcelli ( η / ) [37] was also considered, but the scaling is found to be muchtoo strong.The three fluid theory [19] utilises some ordering assumptions in its derivation; firstly that the electron beta is ofthe order of the mass-ratio ( √ β e ∼ (cid:112) m e /m i ), secondly we are in the low collisionality limit and finally, that all lengthscales are larger then the electron gyro-radius, both of which also hold in our calculations. This model was also derivedwithout a background pressure gradient, however the scaling from our result seems largely invariant to this, as seen inthe further two curves in Fig. 5. The η / scaling found here requires a large ∆ (cid:48) in its derivation. This requirement,in turn, invalidates the constant- ψ approximation. To achieve a constant perturbed magnetic flux across the singularlayer, diffusive processes must act to smooth out any perturbations, thus the growth rate of the mode must be slowerthan the time scale of the diffusive process or the magnetic island must be very narrow [21]. A rough estimate, usingthe resistive growth rate, states that the growth rate τ A γ (cid:28) S − / for the constant- ψ approximation to hold. Thisequality is never really satisfied in the present simulations. In Fig 3 the radial profile of ∂A || /∂ψ shows that throughthe singular layer A || can not be considered constant. Direct calculation of ∆ (cid:48) in this case is difficult in toroidal −3 −2 −1 −2 Collision frequency, ν * [ ν ei q/ ε ] G r o w t h r a t e , γ [ v t h i / R ] −3 −2 −1 −12−10−8−6−4−20 x 10 −3 Collision frequency, ν * [ ν ei q/ ε ] F r equen cy , ω [ v t h i / R ] R/L n = 0R/L n = 1.0 γ collless (R/L n = 0.0) γ = ν ei γ ~ ν ei3/5 γ∼ν ei3/5 γ ~ ν ei1/7 γ collless (R/L n = 1.0)m=2,3 R/L n = 0R/L n = 0 No driftsR/L n = 1.0 ω collless ω e* , diamagnetic freq (R/L n =1.0). ω collless (R/L n = 1.0)m=2,3 FIG. 4. The (left) growth rates and (right) mode frequencies of an m = 2 , n = 1 tearing mode with aspect ratio R/a = 3 in ahydrogen plasma. These have been calculated for purely electron-ion pitch angle collision operator with no background pressuregradient (black circles) and with a background density gradient (red circles). The left hand panel also shows the analyticalresult derived by Drake and Lee for a collisional tearing mode ( γ ∝ ν / ei , green dashes) and semi-collisional mode( γ ∝ ν / ei ,blackdashed line). geometry where the resonant layer width is not precisely defined, as such it can be merely estimated qualitatively.The kinetic calculation of [18] utilises the constant- ψ approximation, as such, with the parameters considered here, ascaling with η / is not expected (Of note, an assumption of ∆ (cid:48) → ∞ is utilised in [20] for an m = 1 perturbation,where the constant- ψ approximation can not be made).The right hand panel of Fig. 4 shows the frequency of the tearing mode as a function of the normalised collisionality.We note that the tearing mode has a small but finite rotation frequency in the absence of a background density ortemperature gradient. Kinetic analysis of the classical tearing mode [17, 18] in slab geometry has shown that, in theabsence of equilibrium pressure gradients, the mode has no frequency ( ω = 0) which is in contrast to our calculatedvalues. Our calculation is a more complete description, performed in three dimensional geometry, and thus includestrapping and curvature drift effects. By artificially switching off the terms related to particle drifts in the gyro kineticequation, it can be seen that this residual frequency is indeed a toroidal effect. When curvature effects are removed,the mode frequency is found to be essentially zero (Right hand panel of Fig 4). At high collisionalities, (i.e ν ∗ > R/L n ) and temperaturegradients ( R/L T ). The same scaling is seen as in Fig. 4, however it is apparent that at high collisionality the modestarts to stabilise, apparent by the reduction of the growth rate.Our result shows that the presence of larger pressure gradients causes a stabilisation of the mode at lower colli-sionalities. This is further evident in Fig. 6 where the growth rate and mode frequency is shown as a function ofthe equilibrium density gradient at fixed collisionality. Initially here the density gradient has a destabilising effect,causing the growth rate to increase. At higher gradients ( R/L n > . −4 −2 −0.04−0.03−0.02−0.0100.01 Normalised collision freq., ν * [ ν ei q/ ε ] M ode , f r equen cy , ω [ v t h i / R ] −4 −2 −2 G r o w t h r a t e , γ [ v t h / R ] η RLn=1 RLt= 2RLn=1, RLt= 0RLn= 1, RLt=1
FIG. 5. The (top) growth rates and (bottom) mode frequencies of an m = 2 , n = 1 tearing mode as a function of thenormalised collision frequency at the rational surface in the presence of a background temperature and density gradients.Horizontal dashed lines (top) represent the collisionless growth rates, while solid lines (bottom) represent the analytic modefrequencies as calculated by ω = ω n ∗ +0 . ω T ∗ in the collisionless limit. Dot-dashed lines are the corresponding analytical valuesfor the semi-collisional frequency, ω = ω ∗ n + 5 / ω ∗ T . stabilising the TM [40]. This effect is a function of the normalised beta, ˆ β = β e / L s /R ) ( R/L n ) where L s = Rq/ ˆ s is the magnetic shear length, and a second parameter, a normalised collisionality, C = 0 . L s /L n ) ( ν ei /ω ∗ e )( m e /m i )which reduces to 0 . δ c /ρ i ) where δ c is the semi-collisional resonant layer width. In the simulations presented it ispossible that δ c ∼ ρ i ([43] assume that ρ i > δ c ).It has been shown analytically [40] that as the parameter ˆ β approaches unity the tearing mode becomes entirelystable. The results presented here have values significantly lower, ( ˆ β ∼ . R/L n ) where ˆ s = 0 .
65 at the q = 2rational surface) but the stabilising effects of pressure gradients still become apparent. Analytical work has shownthat toroidal geometry has a stabilising effect on the tearing mode[41], however, at a low beta, as is used in this work,these effects will be smallThe mode frequency of the tearing mode is plotted in Fig. 5 in the presence of a density and temperature gradientas a function of collisionality along with the corresponding collisionless asymptotes, while in Fig. 6 as a function ofdensity/temperature gradient at fixed collisionality. Kinetic analysis of the mode shows that it rotates at the electrondiamagnetic frequency related to the density and temperature gradient at the resonant layer when they are present.In GKW units the diamagnetic frequency is given by ω ∗ e = − ( R/L n )( k th ρ i ) / − . v thi /R , for R/L n = 1 . R/L T = 0 .
0. When a background temperature gradient is also considered then, in the collisionless limit, the modefrequency is given by the expression, ω = ω ∗ n + 1 / ω ∗ T [18], where ω ∗ T = − / k θ ρ i ) R/L T which is the electrondiamagnetic frequency associated with the temperature gradient, these prediced values are also plotted and theirscalings have been plotted in the right hand panel of Fig. 6. In both cases showing a good agreement. We also see1 −3 G r o w t h r a t e , γ [ v t h i / R ] F r equen cy , ω [ v t h i / R ] Logarithmic density (temperature) gradient, R/L n (R/L T ) ω *n +0.5 ω *T ω *n ν * =0.1 R/L n =1.0 ν * =0.0 R/L T =0 ν * =0.1 R/L T =0 FIG. 6. The growth rate (left) and frequency (right) of the mode as a function of the equilibrium density gradient at theresonant surface (
R/L n , black) or equilibrium temperature gradients ( R/L T , red). The solid grey line is the electron diamagneticfrequency when considering only R/L n while red is the Drake and Lee prediction of the collisionless mode frequency ω ∗ n +0 . ω ∗ T . that the residual frequency as previously mentioned, is significantly smaller than the electron diamagnetic frequencywhen the mode evolves in the presence of pressure gradients producing the so called drift-tearing mode.Plotted in the bottom panel of Fig. 5 are the mode frequencies with (horizontal solid lines) the correspondingdiamagnetic frequency as given by the above expression, giving good agreement in the collisionless limit. It is alsofound that the frequency never reaches the Drake and Lee semi-collisionless expression of, ω = ω ∗ n + 5 / ω ∗ T whichare also plotted as horizontal dot-dashed lines, although a qualitative agreement as a function of the pressure gradientis seen i.e. with a larger background pressure gradient a larger maximal frequency is seen. It is also noted that thecollisionality threshold at which the mode changes the sign of its rotation is a function of the pressure gradient, witha higher threshold seen at higher background gradients. While the mode start to be stabilised at a lower collisionality,the opposite relation. VII. CONCLUSIONS
Using the global version of the gyrokinetic code, GKW, the tearing mode has been studied in the linear regimefor experimentally relevant paramerters in three dimensional, toroidal geometry for the first time. The collisional-ity, equilibrium temperature, density and current gradients were varied and their effects on the mode stability andfrequency were analysed. The main results can be summarised as follows: • The tearing mode growth rate scales with the resistivity close to the η / scaling as calculated by Fitzpatrick[19] in the semi-collisional, high ∆ (cid:48) regime as opposed to the semi-collisional scaling of η / found in the kinetictreatment of Drake and Lee. • The stability of the mode follows closely the stability calculated in cylindrical geometry performed by Wessonet. al. [42]. However, our analysis differs in that multiple rational surfaces may be present within the domainand so allows double tearing modes and toroidal coupling of modes, a process which is known to be destabilising,as such more modes are excited in our global toroidal treatment. • The mode is seen to rotate at the electron diamagnetic frequency when a pressure gradient is present, as ischaracteristic for the drift-tearing mode. However, this is a function of the collisionality and the sign of therotation frequency is seen to change at a sufficienty high collisionality.2 • A residual non-zero, mode frequency, is found for the classical tearing mode in three dimensional toroidalgeometry, even when equilibrium pressure gradients are not included. Suppressing toroidal drift effects recoversa zero mode frequency. • Pressure gradient stabilisation of the mode is apparent, even though β e which determines the effective stabili-sation, is small.Linear mode evolution is valid only when the generated magnetic island is small compared to the resonant layerwidth, which in these simulations were smaller than the ion gyroradius. However, the mode frequency is important forthe non-linear evolution of the magnetic island and in particular the polarisation current drive of the tearing mode,and its neoclassical modification, the Neoclassical tearing mode, particularly when the magnetic island is small [24]. ACKNOWLEDGMENTS
A part of this work was carried out using the HELIOS supercomputer system at Computational Simulation Cen-tre of International Fusion Energy Research Centre (IFERC-CSC), Aomori, Japan, under the Broader Approachcollaboration between Euratom and Japan, implemented by Fusion for Energy and JAEA.The authors would like to thank the Lorentz Center, Leiden for their support. Useful and productive conversationswith C.S. Chang, C. Hegna, G. Huismanns and other attendees of the meeting, “Modelling Kinetic Aspects of GlobalMHD Modes” are gratefully acknowledged. [1] D. Biskamp
Magnetic reconnection in plasmas. 3rd Edition
Cambridge University Press (2005)[2] W. A. Newcomb, Annals of Physics
459 (1963)[4] H. P. Furth, M. N. Rosenbluth, H. Selberg, Phys. Fluids
291 (2009)[6] F. L. Waelbroeck, Nucl. Fusion
899 (1986)[8] C.C. Hegna, Phys. Plasmas
248 (1996)[10] O. Sauter et. al , Phys. Plasmas , 075010 (2009)[16] W.A. Hornsby et.al., Euro. Phys. Lett. , 1778 (1975)[18] J. F. Drake, Y. C. Lee, Phys. Fluids , 1341 (1977)[19] R. Fitzpatrick, Phys. Plasmas , 042101 (2010)[20] A. Y. Aydemir, Phys. Fluids B. , 112106 (2011)[26] B. N. Rogers, S. Kobayashi, P. Ricci, W. Dorland, J. Drake et al. , Phys. Plasmas , 092110 (2007)[27] W. Wan, Y. Chen, S. E. Parker, Phys. Plasmas , 122104 (2012)[29] A.G. Peeters et al (To be submitted)[30] A.G. Peeters, Y. Camenen, F.J. Casson, W.A. Hornsby, A.P. Snodin, D. Strintzi, and G. Szepesi, Comp. Phys. Comm. , 2649 (2009)[31] X. Lapillonne, S. Brunner, T. Dannert, S. Jolliet, A. Marinoni, L. Villard, T. Gorler, F. Jenko, and F. Merz, Phys. Plasmas, , 032308, (2009).[32] S. Hamada, Kakuyugo Kenkyu , 542 (1958)[33] C.F.F. Karney, Comp. Phys. Reports (1986) 183[34] R. Fitzpatrick, P. G. Watson and F. L. Waelbroeck, Phys. Plasmas [36] J. Wesson, Tokamaks (Cambridge University Press, Cambridge, 1987)[37] F. Porcelli, Phys. Fluids
577 (1988)[40] J. F. Drake, T. M. Antonsen Jr., A. B. Hassam and N. T. Gladd, Phys. Fluids
875 (1975)[42] R.J. Hastie, A. Sykes, M. Turner and J.A. Wesson, Nucl. Fusion,54