The linear tearing instability in three dimensional, toroidal gyrokinetic simulations
Abstract
Linear gyrokinetic simulations of the classical tearing mode in threedimensional toroidal geometry were performed using the global gyro kinetic turbulence code, GKW . The results were benchmarked 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, collisionality and the pressure gradients. Both collisionless and semicollisional tearing modes were found with a smooth transition between the two. A residual, finite, rotation frequency of the mode even in the absense of a pressure gradient is observed which is attributed to toroidal finite Larmorradius effects. When a pressure gradient is present at low collisionality, the mode rotates at the expected electron 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 semicollisional regime, closely following the semicollisional scaling found by Fitzpatrick. The stability of the mode closely follows the stability using resistive MHD theory, however a modification due to toroidal coupling and pressure effects is seen.
pacs:
I Introduction
A current density profile in a plasma can lead to instabilities that are able to change the topology of the magnetic field via the process of magnetic reconnection (1). Resistive instabilities, and in particular tearing instabilities (2); (3); (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 a detrimental 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 and therefore the tearing mode, and specifically, the neoclassical tearing mode (NTM) (7); (8) is expected to be the limiting factor on its operational plasma (10).
Much work has been performed numerically in studying tearing mode stability in both cylindrical and toroidal geometry (11), using the Magnetohydrodynamic (MHD) framework. However for large current and future tokamaks a fully kinetic description of the physics in the resonant layer is needed for both the linear and nonlinear regimes.
Linear theory is generally considered to be valid when the generated magnetic island has a width that is smaller than 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 on the evolution of the tearing mode indicates a turbulent modification of the resonant layer, extending the suitability of linear theory to larger island sizes(14), while it has also been seen that a magnetic island can have a profound effect on turbulence and transport (15); (16).
Within the singular layer (Parallel wavevector of the mode, ) the assumptions of ideal MHD break down and magnetic reconnection may take place. For a collisionless mode the mechanism of reconnection is the electron inertia, 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 a classical collisional tearing (3) or, more relevant for present day tokamaks, a semicollisional mode (18); (19). The mode frequency, and in particular, its direction is vital in predicting the nonlinear stabity of the NTM. Specifically the sign of the mode frequency, whether the mode rotates in the ion or electron diamagnetic direction, determines whether the 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 threedimensional toroidal geometry using the gyrokinetic framework of equations. Parameters appropriate to modern and future Tokamak 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 has studied two dimensional magnetic reconnection in a highly magnetised system (27); (25); (26) 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 how the tearing mode is implemented within the gyrokinetic framework. Sec. IV outlines the parameters used, and the resolution requirements. Sec. V and VI presents the scaling of the tearing mode growth rate and mode frequency with collisionality, and with the variation of the background density, temperature and current profiles.
Ii Gyrokinetic framework
The selfconsistent treatment of the tearing mode requires a radial profile in geometry and thermodynamic quantities. We use here the global version of the gyrokinetic code, GKW (29). This is advancement on the already widely used fluxtube 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 assumptions used in driving the tearing instability.
The delta approximation is used. The distribution function is split into a background and a perturbed distribution . The final equation for the perturbed distribution function , for each species, , can be written in the form
(1) 
where is the source term which is determined by the analytically known, background distribution function, is the magnetic moment, is the velocity along the magnetic field, is the magnetic field strength, m and Z are the particle mass and charge number respectively. Here, is used to absorb the time derivative of the parallel vector potential which enters the equations through Ampères law.
The velocities in Eq. (1) are from left to right: the parallel motion along the unperturbed field (), the drift motion due to the inhomogeneous field (), and the motion due to the perturbed electromagnetic field (, where ). The latter is the combination of the velocity () and the parallel motion along the perturbed field line (). Here, the angled brackets denote gyroaveraged quantities.
The background is assumed to be a shifted Maxwellian (), with particle density () and temperature () and toroidal flow velocity ()
(2) 
which determines the source term,
(3)  
(4) 
The thermal velocity , and the major radius () are use to normalise the length and time scales. Using standard gyrokinetic ordering, the length scale of perturbations along the field line () are significantly longer than those perpendicular to the field (). Here, is the normalised ion Larmor radius (where and ).
The electrostatic potential is calculated from the gyrokinetic quasineutrality equation,
(5) 
where the first term represents the perturbed charge density and the second the polarisation density, which is only calculated from the local Maxwellian. The dagger operator represents the conjugate gyroaverage operator. Similarly the parallel vector potential is calculated via Ampères law,
(6) 
The gyroaverage is then calculated as a numerical average over a ring with a fixed radius equal to the Larmor radius using a 32 point interpolation defined by,
(7) 
where is the gyroangle, is the position of the gyrocenter and is the gyroradius. This gyroaverage is used in both the evolution equation of the distribution function, as well as in the quasineutrality and Ampére equations and on both the electron and ion species. The nonlinearity () is neglected in this study and its effects on magnetic island evolution and saturation will be left for a future publication.
GKW uses straight field line Hamada (32) coordinates () where is the coordinate along the magnetic field and is the generalised toroidal angle and is the normalised radial coordinate. For circular concentric surfaces (31), the transformation of poloidal and toroidal angle to these coordinates is given by (30) (. The wave vector of the island is where is the toroidal mode number. GKW uses a Fourier representation in the binormal direction, perpendicular to the magnetic field. The radial and parallel directions are treated using a high order finite difference scheme so as to include global profile effects utilising Dirichelet boundary conditions in the radial direction, and open field line boundaries in the parallel direction. An example of the radial profiles 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 be resistive in nature. In GKW the full linearised Landau collision operator with momentum and energy conserving terms is implemented. However here, unless stated otherwise, we utilise only the pitchangle scattering part. The differential part of the operator has the following form (33)
(8) 
Here the sum is over all species and denotes the particle pitchangle ().
The coefficients can be obtained from the literature (33) and written in normalised form, suppressing the subscript N for the normalised quantities:
(9) 
where and are the standard definition of the error function and its derivative. Finally, the normalised collision frequency, is given by
, is the scattering species number density, Z is the relative charge of the species and is the Coulomb logarithm. and are the velocities of the scattered and scattering species respectively.
Iii Tearing mode implementation
The tearing mode instability is driven by a nonhomogeneous current density. In our implementation this is introduced by applying an electron flow profile in the equilibrium. The electron flow velocity is calculated selfconsistently from the imposed qprofile 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, , i.e. that the electron flow velocity is larger than the ion flow and that it is a flux function. The current gradient is therefore related to the gradient in the electron flow, , this profile is introduced to the code where is written in equation 4.
We assume to be significantly smaller than the electron thermal velocity, thus the only term in Eq. 4 of interest is,
The electron profile is related to the qprofile via Ampéres law, , which in circular geometry and assuming a low aspect ratio, becomes,
(10) 
where the function is approximated to be close to unity and subscript zeros represent constant values. We have used the definition and therefore the poloidal magnetic field, . The gradient of the background current is given by the expression,
(11) 
In this paper we make use of the Wesson analytic form of the current profile (36) which is defined by the expression,
(12) 
which, in turn, gives a profile of the form,
(13) 
where is the safety factor at . This is related to the on the axis, , by , where is an integer that determines that peaking of the current gradient.
The density and temperature profiles have the radial form,
where and are the logarithmic density and temperature gradients at the reference radius, , is the profile width and is the rising width.
Iv Parameters
In this section we study the linear mode structure, a twodimensional example of which is plotted in Fig. 1, showing the radially elognated eigenfunction of the parallel vector potential (left panel) and the electrostatic potential () which is highly localised around the singular layer at the rational surface (dashed line). The growth rates and mode rotation frequencies are calculated for the following parameters unless stated otherwise in the text:
Simulation parameters.  

Aspect ratio, R/a  3 
Electron beta,  
Current peaking parameter,  1 
Mass ratio,  1836 
q at edge,  2.99 
Mode number,  0.0255 
0.002815  
4/8 
Care must be taken in resolving the relevant scale lengths in each of the domain directions. The number of grid points in the poloidal, parallel velocity and directions were, , , respectively which were found to give converged growth rates to within a couple of percent. A single toroidal mode is considered, namely with the mode number, the number of resonant poloidal modes within the domain is determined by the profile (); as such the possibilty of mode coupling effects remain, particularly the toroidal coupling of modes with varying poloidal mode numbers (m) which is known to have a destabilizing effect (38). The radial resolution required is dependent on the physics of the resonant layer. In the collisionless tearing mode the relevant scale length is the electron skin depth, , where is the electron plasma frequency and is the speed of light. The skin depth is, The is the electron defined by , where and are equilibrium density and temperature and is the magnetic field strength at the outboard midplane. For collisional modes, the resistivity defines the width of the layer. Unless otherwise stated, we have used, radial grid points between and gives us three grid points within the skin depth for a normalised gyroradius of, and a of , in the collisionless case without a background pressure gradient. Increasing collisionality and including a background pressure gradient has the effect of increasing the resonant layer width and with it the number of grid points within the layer. A hydrogen mass ratio plasma is studied due to a slightly more relaxed skindepth resolution requirement over a deuterium mass ratio plasma while still maintaining applicability to experimental conditions.
In Fig. 2 an example of the radial eigenfunctions are plotted and compared with the result from an MHD shooting calculation performed in the cylindrical tokamak limit (). The safety factor and magnetic shear profiles used in this calculation are found in the second panel. The shift in the peak of the parallel vector potential with respect to the cylindrical calculation is consistent with previous work performed using an ideal MHD code (11). Plotted is also the radial eigenfunction for a simulation using an aspect ratio of , closer to the cylindrical limit, which shows that the eigenfunction tends toward this cylindrical form as expected. The radial density and geometry profiles used are shown in middle panel. The third panel of Fig. 2 shows the second derivative of the radial profile of the vector potential. This shows clearly the position of the singular layer, the position where ideal MHD breaks down and a discontinuity in the derivative is produced. The bottom panel also shows the radial profile of the electrostatic potential (), showing the highly localised electric field at the rational surface, and the perturbed ion density profile ().
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, . This is analogous to the diagram shown in (42); (36), which shows the stability of the , mode in cyclindrical geometry, denoting the toroidal mode number, while is the poloidal mode number.
The analysis in (36) considered a single , mode, the outlines of which are also plotted in overlay (black lines). However we are unable to decompose the modes in this manner and consider, for example a , mode in isolation. The different coloured points denote the dominant mode within a particular simulation determined by 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 , mode within the domain. However, the results are complicated by the toroidal coupling of modes. Toroidal coupling has the effect of destabilising modes that would otherwise be stable (38); (39), for example in lower panel of Fig. 3 we see the radial eigenstructure of a combined and a mode, the latter of which would otherwise be stable. The point at , has only a 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 mode showing a general increase over the pure mode.
The growth rate of the tearing mode is a function of the well known parameter (3),
(14) 
across the singular layer, whose position is denoted by and upper and lower boundaries by and . This parameter represents the discontinuity in and measures whether mode growth is energetically favourable ( for the mode to grow). As the growth rate of the mode is directly proportional to , measuring the growth rate allows us to study the mode stability. As such the above difference in growth rate between the coupled tearing mode and the pure shows that the toroidal coupling in this case has a destabilising effect (Increasing ), increasing the growth rate by a factor of 2.
It can be seen that points along the line are stable to the tearing mode. For the collisionality used in the figure (, ), a toroidal MHD analysis (Fig.8 in (42)) predicts that points along this line would indeed be stable, and the results from GKW are consistent with this, however, the mode is stabilised at a much lower Lundqvist number, than is predicted in (42). is defined as the ratio of the Alven time () and the resistive time scale , where .
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 the trappingdetrapping rate, ), corresponding to a Ludquist number, S, range between and .
A scan with no background density gradient (Black circles) and with a background gradient scale length of (Red circles) is presented. Firstly, at low collisionality the collisionless tearing mode is observed, their growth rates are plotted as horizontal dashed lines which are the low collisionality asymptotic values. These asymptotic values were checked for convergence by doubling the radial resolution (), the growth rate changed by less than 3% from the value using .
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 and marks the boundary between the collisionless and semicollisional mode regimes, at this point collisions limit the electron response, which allows them to be accelerated, and broadens the resonant layer width. The semicollisional theory only holds in the limit when and the transition between semicollisional and collisionless modes is clearly seen in our simulations. The growth rate obtained deviates significantly from the collisionless rate once the growth rate is greater than the electronion collision frequency.
Plotted are the analytic scalings for the resistive (Green dashes, ) and the semicollisional tearing mode (Black dashes, ) obtained using a kinetic treatment by Drake and Lee (18). However both of these give significantly stronger scalings than our calculated values. Good agreement, particularly when a density gradient is present, is found when we consider the semicollisional scaling at high as given by the three fluid, two dimensional slab calculation by Fitzpatrick (19). This predicts a scaling of the growth rate with (Light grey dashes). This corresponds to regime III in Fig.1 of (19) and was also found in a semicollisional study of modes (20). The viscoresistive scaling as calculated by Porcelli () (37) was also considered, but the scaling is found to be much too strong.
The three fluid theory (19) utilises some ordering assumptions in its derivation; firstly that the electron beta is of the order of the massratio (), secondly we are in the low collisionality limit and finally, that all length scales are larger then the electron gyroradius, both of which also hold in our calculations. This model was also derived without a background pressure gradient, however the scaling from our result seems largely invariant to this, as seen in the further two curves in Fig. 5. The scaling found here requires a large in its derivation. This requirement, in turn, invalidates the constant approximation. To achieve a constant perturbed magnetic flux across the singular layer, diffusive processes must act to smooth out any perturbations, thus the growth rate of the mode must be slower than the time scale of the diffusive process or the magnetic island must be very narrow (21). A rough estimate, using the resistive growth rate, states that the growth rate for the constant approximation to hold. This equality is never really satisfied in the present simulations. In Fig 3 the radial profile of shows that through the singular layer can not be considered constant. Direct calculation of in this case is difficult in toroidal 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, a scaling with is not expected (Of note, an assumption of is utilised in (20) for an 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 or temperature gradient. Kinetic analysis of the classical tearing mode (18); (17) in slab geometry has shown that, in the absence of equilibrium pressure gradients, the mode has no frequency () which is in contrast to our calculated values. Our calculation is a more complete description, performed in three dimensional geometry, and thus includes trapping and curvature drift effects. By artificially switching off the terms related to particle drifts in the gyro kinetic equation, 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 ) the frequency begins to decrease and eventually changes sign, rotating in the ion diamagnetic direction, something not predicted by previous analyses.
Presented in Fig. 5 is a collisionality scan (over the same range of parameters as Fig. 4) showing the growth rate and mode frequency in the presence of three different combinations of background density () and temperature gradients (). The same scaling is seen as in Fig. 4, however it is apparent that at high collisionality the mode starts 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 collisionalities. This is further evident in Fig. 6 where the growth rate and mode frequency is shown as a function of the equilibrium density gradient at fixed collisionality. Initially here the density gradient has a destabilising effect, causing the growth rate to increase. At higher gradients (), the mode begins to stabilize. There are two key parameters in determining the nature of the tearing instability when considering the effects of equilibrium pressure gradients (40); (43), which are noted for having a stabilising effect on the mode as they produce currents outside the resonant layer. The currents shield the singular layer from the magnetic perturbation, preventing reconnection and thus stabilising the TM (40). This effect is a function of the normalised beta, where is the magnetic shear length, and a second parameter, a normalised collisionality, which reduces to where is the semicollisional resonant layer width. In the simulations presented it is possible that ((43) assume that ).
It has been shown analytically (40) that as the parameter approaches unity the tearing mode becomes entirely stable. The results presented here have values significantly lower, ( where at the rational surface) but the stabilising effects of pressure gradients still become apparent. Analytical work has shown that 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 small
The mode frequency of the tearing mode is plotted in Fig. 5 in the presence of a density and temperature gradient as a function of collisionality along with the corresponding collisionless asymptotes, while in Fig. 6 as a function of density/temperature gradient at fixed collisionality. Kinetic analysis of the mode shows that it rotates at the electron diamagnetic 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 , for and . When a background temperature gradient is also considered then, in the collisionless limit, the mode frequency is given by the expression, (18), where which is the electron diamagnetic frequency associated with the temperature gradient, these prediced values are also plotted and their scalings have been plotted in the right hand panel of Fig. 6. In both cases showing a good agreement. We also see that the residual frequency as previously mentioned, is significantly smaller than the electron diamagnetic frequency when the mode evolves in the presence of pressure gradients producing the so called drifttearing mode.
Plotted in the bottom panel of Fig. 5 are the mode frequencies with (horizontal solid lines) the corresponding diamagnetic frequency as given by the above expression, giving good agreement in the collisionless limit. It is also found that the frequency never reaches the Drake and Lee semicollisionless expression of, which are also plotted as horizontal dotdashed lines, although a qualitative agreement as a function of the pressure gradient is seen i.e. with a larger background pressure gradient a larger maximal frequency is seen. It is also noted that the collisionality threshold at which the mode changes the sign of its rotation is a function of the pressure gradient, with a 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 regime for experimentally relevant paramerters in three dimensional, toroidal geometry for the first time. The collisionality, equilibrium temperature, density and current gradients were varied and their effects on the mode stability and frequency 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 semicollisional, high regime as opposed to the semicollisional scaling of found in the kinetic treatment of Drake and Lee.

The stability of the mode follows closely the stability calculated in cylindrical geometry performed by Wesson et. al. (42). However, our analysis differs in that multiple rational surfaces may be present within the domain and 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 is characteristic for the drifttearing mode. However, this is a function of the collisionality and the sign of the rotation frequency is seen to change at a sufficienty high collisionality.

A residual nonzero, mode frequency, is found for the classical tearing mode in three dimensional toroidal geometry, even when equilibrium pressure gradients are not included. Suppressing toroidal drift effects recovers a zero mode frequency.

Pressure gradient stabilisation of the mode is apparent, even though which determines the effective stabilisation, is small.
Linear mode evolution is valid only when the generated magnetic island is small compared to the resonant layer width, which in these simulations were smaller than the ion gyroradius. However, the mode frequency is important for the nonlinear 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).
Acknowledgements.
A part of this work was carried out using the HELIOS supercomputer system at Computational Simulation Centre of International Fusion Energy Research Centre (IFERCCSC), Aomori, Japan, under the Broader Approach collaboration 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 conversations with C.S. Chang, C. Hegna, G. Huismanns and other attendees of the meeting, “Modelling Kinetic Aspects of Global MHD Modes” are gratefully acknowledged.Footnotes
 preprint: Preprint: do not distribute
References
 D. Biskamp Magnetic reconnection in plasmas. 3rd Edition Cambridge University Press (2005)
 W. A. Newcomb, Annals of Physics 10 232267 (1960)
 H. P. Furth, J. Killeen, M. N. Rosenbluth, Phys. Fluids 6 459 (1963)
 H. P. Furth, M. N. Rosenbluth, H. Selberg, Phys. Fluids 16 1054 (1973)
 E.G. Zweibel and M. Yamada, Annu. Rev. Astron. Astrophys. 47 291 (2009)
 F. L. Waelbroeck, Nucl. Fusion 49 104025 (2009)
 R. Carreras, R.D. Hazeltine, M. Kotschenreuther, Phys. Fluids 29 899 (1986)
 C.C. Hegna, Phys. Plasmas 5 1767 (1998)
 H.R. Wilson, J.W. Connor, R.J. Hastie, and C.C. Hegna, Phys. Plasmas 3 248 (1996)
 O. Sauter et. al, Phys. Plasmas 4 1654 (1997)
 Y. Nishimura, J. D. Callen, C. C. Hegna, Phys. Plasmas 5 4292 (1998)
 P. H. Rutherford, Phys. Fluids 16 1903 (1973)
 R.J. Hastie, F. Militello and F. Porcelli, Phys. Rev Lett. 95 065001 (2005)
 W. A. Hornsby et al. (Submitted) arXiv:1403.1520
 E. Poli, A. Bottino and A.G. Peeters, Nucl. Fusion 49, 075010 (2009)
 W.A. Hornsby et.al., Euro. Phys. Lett. 91 45001 (2010)
 R. D. Hazeltine, D. Dobrott, T. S. Wang, Phys Fluids 18, 1778 (1975)
 J. F. Drake, Y. C. Lee, Phys. Fluids 20, 1341 (1977)
 R. Fitzpatrick, Phys. Plasmas 17, 042101 (2010)
 A. Y. Aydemir, Phys. Fluids B. 3 3025 (1991)
 F.L. Waelbroeck, Phys. Rev. Lett. 70 3259 (1993)
 T. S. Hahm and L. Chen, Phys. Fluids 29 1891 (1986)
 R. Fitzpatrick, F. L. Waelbroeck, and F. Militello, Phys. Plasmas, 13 122507 (2006)
 E. Poli, A. Bergmann and A. G. Peeters, Phys. Rev. Lett. 94 205001 (2005)
 R. Numata, W. Dorland, G. G. Howes, N. F. Loureiro, B. N. Rogers, T. Tatuno, Phys. Plasmas 18, 112106 (2011)
 B. N. Rogers, S. Kobayashi, P. Ricci, W. Dorland, J. Drake et al., Phys. Plasmas 14, 092110 (2007)
 W. Wan, Y. Chen, S. E. Parker, Phys. Plasmas 12 012311 (2005)
 A. Mishchenko and A. Zocco, Phys. Plasmas 19, 122104 (2012)
 A.G. Peeters et al (To be submitted)
 A.G. Peeters, Y. Camenen, F.J. Casson, W.A. Hornsby, A.P. Snodin, D. Strintzi, and G. Szepesi, Comp. Phys. Comm. 180, 2649 (2009)
 X. Lapillonne, S. Brunner, T. Dannert, S. Jolliet, A. Marinoni, L. Villard, T. Gorler, F. Jenko, and F. Merz, Phys. Plasmas, 16, 032308, (2009).
 S. Hamada, Kakuyugo Kenkyu 1, 542 (1958)
 C.F.F. Karney, Comp. Phys. Reports 4 (1986) 183
 R. Fitzpatrick, P. G. Watson and F. L. Waelbroeck, Phys. Plasmas 12 082510 (2005)
 I. Katanuma, T. Kamimura, Phys. Fluids 23 2500 (1980)
 J. Wesson, Tokamaks (Cambridge University Press, Cambridge, 1987)
 F. Porcelli, Phys. Fluids 30 6 (1987)
 B. Carreras, H. R. Hicks and D. K. Lee, Phys Fluids 24 1 (1981)
 J.W.Connor, S.C.Cowley, R.J. Hastie, T.C. Hender, A.Hood and T.J.Martin, Phys. Fluids 31 577 (1988)
 J. F. Drake, T. M. Antonsen Jr., A. B. Hassam and N. T. Gladd, Phys. Fluids 26 2509 (1983)
 A.H. Glasser, J.M. Greene and J.L. Johnson, Phys. Fluids 18 875 (1975)
 R.J. Hastie, A. Sykes, M. Turner and J.A. Wesson, Nucl. Fusion, 17 3 (1977)
 J.W. Connor, R.J. Hastie and A. Zocco, Plasma. Phys. Contol. Fusion 54 035003 (2012)