Drift kinetic theory of neoclassical tearing modes in a low collisionality tokamak plasma: magnetic island threshold physics

A new drift kinetic theory for the plasma response to the neoclassical tearing mode (NTM) magnetic perturbation is presented. Small magnetic islands of width, w≪a (a is the tokamak minor radius) are assumed, retaining the limit w ∼ ρ bi (ρ bi is the ion banana orbit width) to include finite orbit width effects. When collisions are small, the ions/electrons follow streamlines in phase space; for passing particles, these lie in surfaces that reproduce the magnetic island structure but have a radial shift by an amount, proportional to ρϑi/e , where ρϑi/e is the ion/electron poloidal Larmor radius. This shift is associated with the curvature and ∇B drifts and is found to be in opposite directions for V∥≶0 , where V∥ is the component of velocity parallel to the magnetic field. The particle distribution function is then found to be flattened across these shifted or drift islands rather than the magnetic island. This results in the pressure gradient being sustained across the magnetic island for w∼ρϑi and hence reduces the neoclassical drive for NTMs when w is small. This provides a physics basis for the NTM threshold, which is quantified. In Imada et al (2019 Nucl. Fusion 59 046016, and references therein), a 4D drift kinetic non-linear code has been applied to describe these modes. In the present paper, the drift island formalism is employed. Valid at low collisionality, it allows a dimensionality reduction to a 3D problem, simplifying the numerical task and efficiently resolving the collisional boundary layer across the trapped-passing boundary. An improved model is adopted for the magnetic drift frequency. This decreases the NTM threshold, compared to the results shown in Imada et al (2019 Nucl. Fusion 59 046016, and references therein), making it in quantitative agreement with experimental observations, with wc=0.45ρϑi , where w c is the threshold magnetic island half-width, or 2.85ρ bi for the full threshold island width, predicted for our equilibrium.


Introduction
Neoclassical tearing modes (NTM) [1,2] are classified as large scale resistive magnetohydrodynamic (MHD) plasma instabilities [3]. They limit plasma pressure, reducing fusion gain, and occur in the standard ELMy H-mode, anticipated for the ITER baseline scenario, as well as advanced tokamak scenarios. They can arise in tokamak plasmas due to a filamentation of the plasma current density parallel to the magnetic field lines. This filamentation changes the topology of the magnetic flux surfaces, tearing them apart, to form magnetic islands. Experimentally, NTMs are found to be driven unstable when two criteria are satisfied simultaneously: a threshold in poloidal beta and a threshold island width have to be exceeded. Thus, NTMs are usually triggered above a threshold beta by another MHD perturbation (e.g. sawteeth, fishbones, ELMs, error field locked modes, etc) that creates a primary seed island of sufficient amplitude for NTM growth 3 .
According to the conventional theory [1,2], which requires islands to be sufficiently large, the plasma pressure is quickly equalised around the island due to the large particle and heat transport along magnetic field lines. Thus, in the absence of any heat/particle sources from within the island, the plasma pressure gradient inside the island and hence the total plasma pressure in the core are reduced. This flattening of the pressure profile across the island, in turn, leads to a hole in the bootstrap current near the island O-point, providing the filamentation to drive island growth. As the bootstrap current density rises with beta, the island width also grows with beta, resulting in a degradation of confinement [4][5][6][7]. Along with this soft beta limit, NTMs with lower poloidal mode numbers can also lead to plasma disruptions through mode locking. There are a number of NTM control techniques [6]; one of them is to generate microwaves at the electron cyclotron frequency to drive current inside the island and replace the missing bootstrap current. This O-point electron cyclotron current drive (ECCD) scheme has demonstrated complete NTM stabilisation on a number of machines and is to be applied to drive the island width down to mitigate the confinement degradation and/or suppress the NTM in fusion devices such as ITER [6]. However, an issue here is to determine how much of the ECCD current is required for the NTM stabilisation and how localised should it be. This leads to a requirement for a more detailed understanding of the threshold physics.
The NTM magnetic islands can either grow or shrink, depending on the filamentary current density perturbation localised to the rational surface, and parallel to the magnetic field, J ∥ . According to the modified Rutherford 3 There is also a special class of the so called triggerless NTMs that can grow to the saturated state without any preceding MHD activity. theory [6][7][8][9], the time evolution of the island width is obtained from Ampère's law, which can be expressed in the form: where τ R ∼ µ 0 a 2 /η is the resistive diffusion time, η is the local plasma resistivity, w is the island half-width and r s is the radius of the rational surface, where the safety factor q(r s ) ≡ q s = m/n with m/n being the poloidal/toroidal mode number of the NTM. ∆ ′ is the classical tearing mode stability parameter [10][11][12][13]. It arises due to the free energy in the equilibrium current density, and is calculated from ideal MHD as the discontinuity in the perturbed magnetic flux gradient across the rational surface. In Rutherford's original work [13], only the induced current associated with island growth contributes to J ∥ , and this has been explicitly written through the term on the left hand side of equation (1). Adding tokamak neoclassical effects and particle drift physics, denoted by the second term on the right hand side of equation (1) , leads to the so-called modified Rutherford equation (MRE) [7,8]. The integral over q, where q is a triplet of spatial coordinates, represents an integral across the rational surface of the component of J ∥ in phase with the magnetic island. The MRE main contributions come from the bootstrap [1,2], curvature [14,15] and polarisation [16][17][18] currents and are denoted by ∆ bs , ∆ cur and ∆ pol , respectively. The perturbed bootstrap current exists in the banana collisionality regime in a tokamak and is written through a linear combination of the electron/ion density and temperature gradients. For magnetic islands well above the threshold with no heat/particle sinks/sources, the pressure gradient and the bootstrap current perturbation tends to be removed from within the island, and ∆ bs ∼ ε 1/2 (L q /L p ) (β ϑ /w) [9,18] and hence is destabilising 4 . Here ε = r s /R 0 is the inverse aspect ratio with R 0 being the tokamak major radius, β ϑ is poloidal beta; the safety factor and pressure length scales are L −1 q,p = ±∇ r ln q, p > 0 for positive shear. However, the physics is much more complicated for smaller w, relevant for the threshold. According to experimental observations [19,20], small magnetic islands heal themselves. This fact proves the existence of additional physics to provide the tearing mode threshold. One such threshold mechanism originates from the effects of finite radial diffusion [21,22,30] and another arises from finite orbit widths [16,18,[23][24][25][26][27][28][29][30].
Finite orbit width effects are associated with small magnetic islands of width comparable to the trapped ion banana orbit width, when the polarisation current generally plays a role alongside the bootstrap current. The polarisation current is sensitive to the physics of a layer around the island separatrix.
The model we discuss below is not particularly rigorous for this. So we restrict consideration to the situation where the equilibrium radial electric field is zero in the island rest frame, in which case the polarisation current is not expected to be important. The dominant physics is then associated with the bootstrap current and curvature.
∆ cur describes the stabilising curvature contribution as an extension of the linear theory, introduced by Glasser, Greene and Johnson [14,15]. Being an O(ε 2 ) effect, the curvature contribution is generally weak in the large aspect ratio tokamak geometry we consider here. However, in spherical tokamaks the curvature contribution and the bootstrap drive can be comparable [31].
Existing theories of NTMs typically require the island width to be much larger than the ion banana orbit width. There is no analytic theory developed for the neoclassical contributions for w ≲ ρ bi . The MRE form shown in [8] is continued heuristically to a region where w ∼ ρ bi , but there is no rigorous theoretical justification for it. However, in [32] it has been shown that the marginal island width below which the NTM is removed, i.e. dw/dt < 0, is about 2ρ bi in both ECCD and beta rampdown discharges and is about 3ρ bi in [33]. This is exactly the region where the existing theory breaks down. Thus, a new theory is required to determine all the MRE neoclassical contributions allowing the limit of w ∼ ρ bi , which is crucial in providing the NTM threshold island width scaling for ITER and other future tokamak devices.
We will consider the low collision frequency limit, ν ii /ε < k ∥ V Ti and ν e /ε < k ∥ V Te , where ν ii /ν e is the ion/electron collision frequency and V Te/i is the electron/ion thermal velocity. k ∥ = k ϑ w/L s , where L s = Rq/s is the shear length scale with s = (r/q)dq/dr being the magnetic shear, k ϑ = m/r s . We extend our previous results [34][35][36] to treat the electrons with the same drift kinetic formalism that we use for the ions, while considering magnetic islands at rest in the plasma E × B frame (ω = 0, where ω is the island propagation frequency in that frame). Earlier work on the ion bootstrap flow considered small islands of w ∼ ρ bi by solving the drift kinetic equation through a Monte Carlo computational approach [37,38]. This simulation confirmed that the ion density gradient is not removed from the region inside small islands. However, [37,38] focused on the ion response only, omitting the electron response due to the narrowness of the electron banana orbit, and hence neglected the effects of the electrostatic potential that arises from plasma quasi-neutrality. The analytic reduction described here explains the physical origin of the density gradient across the island and provides a new NTM threshold model that arises from both, ion and electron plasma responses. It also provides the self-consistent electrostatic potential, Φ, required to ensure quasi-neutrality. When w ≫ ρ ϑi , the electron and ion distribution functions reproduce the results of the original theory [18] and Φ does not play a major role when the island is at rest in the E × B reference frame. However, when ρ ϑe ≪ w ∼ ρ ϑi , we find that the electron and ion solutions differ significantly near the island, highlighting the importance of deriving Φ self-consistently from plasma quasi-neutrality. Once the plasma responses are found, we proceed to the NTM threshold width calculation determining the total perturbed current density along the field lines.
The remainder of the paper is organised as follows. Section 2 introduces the magnetic geometry and the mode dispersion relation. In section 3 we derive the equations describing the plasma response to the NTM magnetic perturbation. The self-consistent electrostatic potential is found in section 4. The drift magnetic island concept is described in section 5. In section 6 we calculate the neoclassical contributions to the MRE and determine the threshold magnetic island width. A conclusion follows. More detailed information, including benchmarking against other models, can be found in appendices A-C.

Magnetic topology and NTM dispersion relation
A small inverse aspect ratio, circular poloidal cross section tokamak approximation is considered. A triplet of spatial variables {ψ, φ, ϑ} forms an orthogonal set of coordinates 5 , satisfying ∇φ × ∇ψ = rB ϑ ∇ϑ, where ψ is the poloidal flux function, while φ and ϑ are the toroidal and poloidal angles, respectively; r is the tokamak minor radius and B ϑ is the poloidal component of the magnetic field. The equilibrium magnetic field is given by: where I = RB φ depends on the poloidal current, R is the tokamak major radius and B φ is the toroidal component of the magnetic field. As ε ≪ 1 and We employ a low beta approximation and to ensure zero divergence of the total magnetic field, we take a magnetic field perturbation associated with the tearing mode to be of the form: A ∥ is the parallel component of the vector potential connected to the NTM poloidal flux perturbation, δψ, via: with δψ =ψ cos nξ provided a single isolated NTM island is considered. Here ξ is a helical angle in the island rest frame defined as: is the NTM perturbation amplitude with w ψ being the island half-width in ψ space related to w in r space via w = w ψ / (RB ϑ ). q ′ s denotes dq/dψ evaluated at the 5 ∇ψ · ∇ϑ = ϑ ′ R 2 B 2 ϑ with ϑ ′ = ∂ϑ/∂ψ| χ will provide an O(ε 2 ) correction. Here χ is the poloidal variable such that ∇ψ · ∇χ = 0 and {ψ, φ, χ} forms an orthogonal set of coordinates. resonant surface, ψ = ψ s . For further analysis, it is convenient to switch from {ψ, φ, ϑ} to {ψ, ξ, ϑ}. To describe the magnetic island geometry, we introduce a perturbed flux surface function Ω that satisfies B · ∇Ω = 0: with Ω = 1/ − 1 corresponding to the separatrix/O-point of the magnetic island, respectively. To derive the NTM dispersion relation, we address Ampère's law written along the field lines. Projecting out the J ∥ components that are in and out of phase with the magnetic island and integrating Ampère's law across the island yield: Here we have assumed a stationary magnetic island in the plasma E × B reference frame to simplify the analysis below 6 . J ∥ is the ϑ-average of J ∥ . Equations (7) and (8) provide a system to be solved for the threshold magnetic island half-width, w c , and the island propagation frequency, ω, once the perturbed current localised about the island, J ∥ , is obtained. This is to be calculated from the ion and electron distribution functions, which we find in the following sections.

Plasma response
The ion/electron response to the NTM magnetic perturbation is described by the drift kinetic equation that reads: in the island rest frame for each particle species, j. A system of two particle species is addressed: plasma electrons and ions. f j here is to be understood as the gyro-angle independent, leading order distribution function in an expansion in ρ cj /L ≪ 1, where ρ cj is the Larmor radius of species j and L is the characteristic size of the system. All spatial derivatives in equation (9) are to be calculated at fixed magnetic moment, µ = V 2 ⊥ /2B, and kinetic energy, K = V 2 /2, where V is the particle speed. (dµ/dt)∂f j /∂µ is omitted as a higher order correction since dµ/dt = O(ρ cj β/L) 7 . Here ∥ denotes a vector component along the magnetic field lines, and magnetic (∇B and curvature) drift contributions, respectively. ω cj = eZ j B/m j is the cyclotron frequency; eZ j and m j are the particle charge and mass. Φ is the electrostatic potential to be determined from the quasineutrality requirement. C j here is a model integro-differential collision operator chosen as in [18] and reproduced in equation (B2).
Expanding about a Maxwell-Boltzmann equilibrium plasma, we write being the Maxwell-Boltzmann distribution of a species j. The plasma density, n eqm , is related to n 0 by n eqm = n 0 (1 − eZ j Φ/T j ), provided eZ j Φ ≪ T j , where Φ is the electrostatic potential in the island rest frame. V Tj = (2T j /m j ) 1/2 is the thermal velocity of a species. g j describes the perturbation in the particle distribution due to the tearing mode occurrence and also includes the equilibrium neoclassical physics. Seeking a solution localised to the rational surface, we Taylor expand the equilibrium around ψ = ψ s , i.e. where and Φ = Φ ′ eqm ψ=ψs (ψ − ψ s ) + δΦ (prime denotes the derivative with respect to ψ, unless otherwise stated), and thus Φ (ψ s ) = δΦ. Φ eqm is the equilibrium potential in the absence of the island, and δΦ is the perturbation associated with a difference in the electron and ion responses to the magnetic island. For simplicity, below we neglect the equilibrium radial electric field, i.e. Φ ′ eqm = 0, in the threshold calculation. This is equivalent to considering an island which is at rest in the E × B rest frame. The perturbed distribution, g j , then must be linear in ψ far from the island to match to the Maxwellian equilibrium, ∂g j /∂ψ| ψ→±∞ = ∂ ψ f M j (ψ s ), and the electrostatic potential satisfies ∂Φ/∂ψ| ψ→±∞ = 0.
To solve equation (9) for g j , we define a small parameter ∆ = w/a ≪ 1 with the following orderings: eZ j Φ/T j ∼ ∆, δΦ/Φ eqm ∼ ∆, B 1 ∼ ε∆ 2 B 0 and g j /f M j ∼ ∆. Considering equation (27) for electrons and equation (39) for ions from [18], we notice that the dimension of the problem can be reduced by switching from {ψ, ξ, ϑ, µ, V} to {p φ , ξ, ϑ, µ, V}, where p φ = ψ − ψ s − IV ∥ /ω cj is the toroidal component of the canonical angular momentum. IV ∥ /ω cj is the excursion of a particle orbit from the reference flux surface. As w ≪ a, the plasma is toroidally symmetric to leading order and thus we expect the toroidal component of p φ to be approximately constant on a particle orbit. To rigorously demonstrate this, we employ an expansion in ∆ and write g j = α g  j , is independent of ϑ at fixed p φ , i.e. g j (p φ , ξ, µ, V). As we consider the banana collisionality regime, the collision operator on the right hand side of equation (9) has been assumed to be order ∆ smaller than the free streaming.
Proceeding to next order in ∆, we obtain in a low beta limit: B 1 · ∇ϑ and B 1 · ∇ξ have been neglected as higher order terms in the limit of small magnetic islands 8 . We note that ∂/∂ψ, ∂/∂p φ ∼ 1/RB ϑ r when acting on equilibrium quantities and ∂/∂ψ, ∂/∂p φ ∼ 1/RB ϑ w on perturbed quantities (∂p φ /∂ψ = 1 to leading order in ρ ϑj /a). To solve equation (11) for g (0) j , we have to eliminate the term in g (1) j , integrating the equation over ϑ at fixed p φ . This is equivalent to an orbit-averaging procedure. For passing particles, g j is periodic in ϑ and thus we simply integrate over a period in ϑ, assuming g j (−π) = g j (π). Trapped particles oscillate between bounce points, ±ϑ b , defined from λB 0 (ϑ b ) = 1, where V ∥ tends to zero. Thus, for them we integrate between ±ϑ b and sum over σ, where σ = V ∥ / V ∥ and λ = 2µ/V 2 is the pitch angle. As continuity is required at each bounce point, this annihilates the ∂g To employ the collision operator from [18], we switch from {µ, V} to {λ, V; σ} in velocity space. Thus, the velocity space integral and V ∥ become: The trapped-passing boundary in pitch angle space is then located at the inverse of the maximum value of the magnetic field, i.e. λ c = 1/B 0 (1 + ε) for the equilibrium given above. λ ∈ [0, λ c ] for passing and λ ∈ (λ c , λ fin ] with λ fin = 1/B 0 (1 − ε) for trapped particles. Therefore, an orbit-averaged form of equation (11) reads: where are the magnetic and E × B drift frequencies in ξ and radial directions, respectively, to O(ε∆g j ) 9 . Θ denotes the Heaviside step function. The ϑ-averaging operator at fixed p φ is defined as: Φ has been assumed to be periodic in ϑ. Using equation (3), we find A detailed derivation can be found in [39]. 9 See appendix A for more detail.

Plasma quasi-neutrality and electrostatic potential
We adopt a Maxwell-Boltzmann equilibrium and so we obtain:n for the ion/electron density integrating equation (10) over velocity space. Heren j = n j /n 0 , δn j = δn j /n 0 and δΦ = eδΦ/T j provided Z i = 1 and T e = T i (this assumption is maintained throughout the paper unless otherwise stated). δn j is the perturbed density associated with g j . Thus, balancing the electron and ion densities, we find: for the perturbed potential. We adopt an iteration process such that at each iteration step, k, we solve i,e (δΦ (k) ) are ion and electron densities for δΦ = δΦ (k) . The iteration proceeds until δΦ (k+1) = δΦ (k) to some defined accuracy. Thus, in contrast to [37], we allow the restoration of the density/temperature gradient across the magnetic island to be influenced by both ions and electrons.

Drift islands
To analyse equation (14) further, we introduce the following dimensionless system: (note: λ is kept non-normalised and we adopt the convention ψ = 0 at the magnetic axis). Employing the large aspect ratio, circular cross section tokamak approximation, assuming that the fastest radial variation is in perturbed quantities and taking ∂p φ /∂ψ = 1 to leading order in ρ ϑj /a, we derive: from equation (14), where the stream function, S, is Hereω D andĈ j g (0) j denote the normalised forms of ω D and the right hand side of equation (14) [39]. Their explicit representations can be found in appendix B, where we show thatĈ j contains the pitch angle scattering and diffusion in S space. S is to be treated as a new radial coordinate. We note that, in contrast to ψ or p φ , S has a different representation for passing and trapped particles; this will provide a complication for matching at the trapped-passing boundary, λ = λ c .
Equations (20) and (21) complete the transformation from {p φ , ξ, λ, V; σ} to {S, ξ, λ, V; σ}, and the particle distribution function is now to be found as g According to its definition, S is a function of p φ , ξ, λ and V for each σ, and depends on the form of the electrostatic potential, which is, in turn, a function of ψ, ξ and ϑ. For passing particles in the absence of the electrostatic potential, i.e. when the E × B drift effects are ignored, the contours of constant S reproduce the magnetic island structure given by the constant Ω contours (see equation (6)) but have a radial shift by the amount ω DρϑjLq /ŵ +ρ ϑjV∥ , which is proportional to the poloidal Larmor radius. This shift arises from ∇B and curvature tokamak drifts, and asω D is σ-dependent in the passing branch, the shift is in opposite directions for V ∥ ≷ 0. We refer to these island structures in the contours of constant S as drift islands. Inclusion of Φ, in principle, might modify the structure of constant S contours. However, the electrostatic potential calculated iteratively to ensure plasma quasi-neutrality does not add any significant qualitative modifications to the form of S; specifically, the surfaces of constant S in the (x, ξ) plane are island-like for passing and open for trapped particles even with Φ included. Examples for passing particles are shown in figure 1 for (a) ρ ϑi /w = 0.05 ≪ 1 and (b) ρ ϑi /w = 0.4. A similar drift island structure in view of plasma tokamak transport has been identified by Kadomtsev in [40], where the chains of islands much smaller than ρ ϑi but larger than ρ ϑe are considered.
To solve equation (20) for g (0) j as a function of S, we employ weak collisional dissipation. In the reference frame in which the equilibrium radial electric field is zero, this is equivalent to imposing δ i ≡ ν ii /εk ∥ V Ti ≪ 1 for ions and δ e ≡ ν e /εk ∥ V Te ≪ 1 for electrons. Employing the perturbation theory again, and applying an expansion in δ j , we obtain ∂g (0,0) j /∂ξ S,ϑ = 0 to leading order. From this we learn that . Thus, we find that the particle distribution function is constant along these S streamlines in the absence of collisions. Proceeding to next order in δ j and introducing collisions, we derive: where A denotes the coefficient in front of ∂g on the left hand side of equation (20). Equation (22) allows us to reconstruct the actual form of the ion/electron distribution, i.e. its S and λ dependence. To eliminate the term in g (0,1) j , we divide both sides of equation (22) by A and introduce the annihilation operator similar to equation (16) to provide ξ-averaging at fixed S. For open contours (i.e. for passing particles outside the drift island and for trapped particles), we integrate equation (22) over a period in ξ, imposing periodicity at ξ = ±π. For closed contours, i.e. passing particles inside the drift island, we integrate between the ξ-bounce points given by  . The σ = ±1 drift islands are centred around pφ(σ = ±1), denoted by grey vertical lines. The magnetic island is located between them;pφ = ±1 corresponds to the separatrix of the magnetic island.
for each ξ, λ, V and σ, and also sum over the two sides of the closed contours, σ pφ = ±1, where σ pφ is the sign of p φ − p φ0 . Continuity at each ξ-bounce point then eliminates g (0,1) j . Thus, equation (22) reduces to: with the ξ-averaging operator at fixed S being defined as: for passing (S c denotes the drift island separatrix and is to be updated at the end of each iteration in Φ for passing particles) and: for trapped particles. The explicit form of the Ĉ j /A S ξ operator is discussed in appendix B. From the S island formalism we learn that while collisions are neglected in equation (20), the combined effect of the parallel streaming, ∇B and curvature drifts would force the passing particle distribution to be flattened inside these drift islands rather than the actual magnetic island. Introducing collisions at next order fully determines the perturbed distribution function via equation (23) and the results of appendix B. This provides a physics basis for the NTM threshold by reducing the neoclassical drive for NTMs of width w ∼ ρ ϑi when the radial shift in the contours of constant S, compared to Ω, is significant.
A schematic ion distribution function vs.p φ is shown in figures 2(a) and (b) at small and large ρ ϑi /w. At small ρ ϑi /w and hence small radial shift in S, equation (21), a sum of the ion distribution functions over σ = ±1 is found to be flattened inside the magnetic island in the vicinity ofp φ = 0. From equation (12), this would result in flattening of the ion density profile around the magnetic island O-point for ρ ϑi /w ≪ 1. In contrast, when the radial shift of the drift islands compared to the magnetic island becomes significant, the flattening of σ g (0),σ i and hence the density flattening are removed from inside the magnetic island, e.g. ρ ϑi = 7.0 × 10 −3 r s is sufficient to partially restore the density gradient across the magnetic island of width w = 0.02r s . If ρ ϑi /w ≳ 1, the profile will be further steepened across the O-point. This results in the density gradient being restored in the vicinity of the island O-point. The corresponding density profiles are shown in figure 3. We highlight that the gradient inside the magnetic island is a consequence of the drift island structures, and thus is a property of the passing (and not trapped) particles.
For electrons, the radial shift in equation (21) is small as ρ ϑe ≪ ρ ϑi . Hence, the drift island effect is less significant for the electron distribution function. This creates a significant difference in the electron and ion density profiles especially at large ρ ϑi in the absence of the electrostatic potential. Indeed, when ρ ϑi /w ≪ 1, the ion and electron density gradients are both removed from inside the magnetic island due to strong parallel streaming along perturbed flux surfaces, so there the role of Φ is not crucial. In contrast, when ρ ϑi and w are comparable, a non-zero, finite ion density gradient is sustained around the magnetic island O-point, while the electron density gradient is still removed in the absence of any potential due to the strong electron parallel streaming and ρ ϑe ≪ w. However, to maintain plasma quasi-neutrality, the electrostatic potential is required, which adjusts to provide n i ≈ n e . Hence, the ion density steepening at large ρ ϑi is explained by the radial shift in S given by equation (21), while the sustainability of the electron density gradient is associated with the self-consistent electrostatic potential.
The particle flow along the field lines is u ∥j = n −1 0´g j V ∥ dV, and thus σ σg (0),σ i represents the parallel flow moment due to equations (12) and (13). The main contribution to the flow is provided by passing particles due to the summation over σ in the ϑ-averaging operator introduced for trapped particles, equation (16). However, the trapped branch also contributes as the integration here is imposed at fixed ψ, and the trapped particle distribution function is g In figure 4 we show contours of the parallel component of the ion flow, u ∥i , at different ρ ϑi /w in the (ψ, ξ) plane. In  to σ σg (0),σ i inside the island and thus there is some flow that penetrates into the island. This contribution grows with growing ρ ϑi /w and spreads further into the island, which provides the basis for an NTM threshold. The corresponding ϑ-averaged parallel component of the current density that contributes to equation (7) is shown in figure 5 for small and large ρ ϑi /w. The corresponding electrostatic potential is shown in figure 6.

Dissipation layer in the vicinity of the trapped-passing boundary
The perturbative approach we applied above in section 5 to eliminate the ξ-dependence at fixed S breaks down in a dissipation layer, i.e. a narrow region in pitch angle in the vicinity of the trapped-passing boundary, λ = λ c (see figure 7). Here collisional dissipation becomes comparable to parallel streaming, ν ii/e ∂ 2 /∂λ 2 ∼ A ∂/∂ξ| S , due to a steep gradient in λ, and thus a full solution of equations (14)/ (20) is required in the layer, close to λ = λ c . This steep gradient region arises because the different orbit-averaging procedures for passing and trapped particles produce a discontinuity at λ = λ c , and the layer resolves this discontinuity.
Following [18], we impose the matching conditions given by equations (106)-(108) of [18] at the trappedpassing boundary to provide continuity of the particle distribution function and its first λ derivative across the boundary. We note that originally matching is imposed at fixed ψ. However, switching from ψ to S and solving equation (23) at the 0th iteration in Φ, we find g . Inset: a full solution of equation (20) in a collisional layer around λc. The trapped branch solution is σ-independent at fixedpφ due to the summation over σ in equation (16).
p φ /ψ without introducing the layer explicitly. The introduction of this layer allows the particle distribution to vary on S contours and hence provides the freedom to impose the matching conditions, equations (106)-(108) of [18], at λ = λ c .
Taking into account the narrowness of the dissipation layer, we fix all the coefficients in equation (14) at λ p/t ≡ λ c ∓ ϵ, where ε is the width of the layer. This is then equivalent to equation (20) with coefficients evaluated at λ p/t and S being replaced with S =Ŝ + ∂ λ S λ p/t λ − λ p/t , whereŜ = S(λ p/t ). Exploiting the thinness of the layer again, we write: and thus within the layer our equation for g (0) j becomes: Hereν j is to be understood as ν ii /V Ti for ions and (ν ee + ν ei )/V Te for electrons, and a is defined as . Similar to [41], we replace ξ with an angle variable, γ ±/t , defined below, and hence equation (26) can be reduced to a simple diffusion equation: where D ±/t = (2ν j /V)a λ p/t for passing, σ = ±1, and trapped branches. Here we have introduced: for passing particles outside theŜ island and for trapped particles. For passing particles inside theŜ island, γ ±/t has the same features as the angle variable introduced in [41].
angle variable. λ = 0 defines the trapped-passing boundary; λ ≶ 0 corresponds to the passing/trapped region, respectively. In contrast to [18], our layer solution includes both regions inside and outside the magnetic island. Equation (27) allows an analytic solution of the Fourier form: The Fourier coefficients, a ±/t n , b ±/t n (n ≥ 0), are unknown and to be found from matching at λ = 0: H ±/t represents a sum of the drive term/contribution from outside the layer and the 0th harmonic, a ±/t 0 [39]. Equation (31) is a set of three equations for 6N + 3 unknowns, n ∈ [0, N], where N represents the Fourier harmonics retained. Due to a difference in γ ±/t , matching at fixed ψ/p φ cannot be provided in n space in the presence of Φ. However, γ ±/t and n are conjugated variables, and γ ±/t is connected with ξ via equations (28)- (30). Thus, taking a number of points in ξ space N ξ = 2 N + 1 and treating γ ±/t = γ ±/t Ŝ , ξ, V = γ ±/t Ŝ (p φ , ξ, V) , ξ, V , we solve equation (31) numerically for a ±/t n , b ±/t n , providing matching at fixed p φ and ξ. Substituting the obtained Fourier coefficients into the above analytic solution in the layer provides matching across the trappedpassing boundary and is then to be used to provide the boundary condition required to solve for g (0,0) j from equation (23). In previous sections we have outlined the derivation of the orbit-averaged drift kinetic equation to leading order in ∆ for ions and electrons that takes into account the electrostatic potential consistent with the plasma quasi-neutrality condition. Two numerical codes have been developed. One of them finds a solution of equation (14), which is a function of p φ , ξ and λ and keeps collisions at leading order for a full range of λ variation. This is to be referred to as the DK-NTM solution and is explained in [36] and [42]. To introduce the drift islands explicitly for both electrons and ions and to efficiently resolve the collisional dissipation layer, the RDK-NTM solver has been developed to describe low collision frequency regimes [39]. It finds a solution of the reduced orbit-averaged drift kinetic equation to leading order in ∆, i.e. equation (14) in the dissipative layer and equation (23) outside the layer with the electrostatic potential consistent with plasma quasi-neutrality. The RDK-NTM numerical algorithm is described in [39]. The numerical results presented in this paper are provided by the RDK-NTM code. In appendix C we compare results from both models.

Neoclassical drive for NTMs
We can now return to equation (1) and consider how J ∥ contributes to the evolution of the magnetic island width. We note that equation (7) is equivalent to equation (1) provided a single isolated stationary NTM magnetic island is considered. Thus, for the stationary island, the classical tearing mode stability parameter, ∆ ′ , is balanced against the sum of all the neoclassical contributions, ∆ ′ + ∆ neo = 0, where: Using the obtained ion/electron distribution function, we calculate J ∥ = j eZ j u ∥j . Defining the polarisation current density as the part of the parallel current density perturbation that flux surface averages to zero, we write: for the sum of the bootstrap and curvature contributions and hence, for the polarisation term. Here the ξ-averaging operator at fixed Ω is defined as: similar to equation (24). As we mentioned earlier, we focus on a large aspect ratio, circular cross section tokamak approximation here, and employ our 3D RDK-NTM code to analyse the threshold physics. Since ∆ cur = O(ε 2 ), it does not provide a significant contribution to the threshold obtained in this work. In particular, ∆ bs + ∆ cur reduces to the bootstrap current contribution for magnetic islands in the limit of large widths, w ≫ ρ ϑi .
In figure 8 we plot ∆ bs + ∆ cur against w. In the limit of w ≫ ρ ϑi , ∆ bs + ∆ cur is inversely proportional to w, which is expected from the existing analytic theory, e.g. equation (85) of [18]. When w tends to zero, ∆ bs + ∆ cur becomes negative, providing a threshold for NTMs, i.e. a value of w below which the mode is stable, ∆ bs + ∆ cur < 0. This value is denoted by w c and is to be referred to as the critical magnetic island half-width. w c is different for each ion poloidal Larmor radius and hence can be defined as a function of ρ ϑi . This kind of behaviour at w ∼ ρ ϑi is the direct result of the inclusion of the drift islands in our model and is in agreement (qualitative and quantitative) with experimentally observed self-healing of small magnetic islands below the threshold (e.g. [33]). A straight line fit shown in the inset of figure 8. provides w c = 0.46ρ ϑi and in terms of the ion banana width w c = 1.47ρ bi at ε = 0.1 (here w c , ρ ϑi and ρ bi are normalised to r s ). w c /ρ ϑi as a function of ε is shown in figure 9, which indicates that the banana width is the appropriate measure. Recall that experimentally it was found that the full critical magnetic island width is between (2, 3)ρ bi , so this result is in the right vicinity. In previous work [36] we found a much larger threshold w c ≃ 9ρ bi , but there we had a larger magnetic drift frequency-the value chosen here, with L −1 B = 0, is closer to the experimental case 10 . This highlights an additional sensitivity of w c to L B .
When the island propagation frequency is close to zero, the polarisation current contribution is expected to be small. Indeed, we find that inclusion of ∆ pol at ω = 0 does not shift the threshold significantly falling slightly to w c = 1.41ρ bi . The latter is obtained in the conventional tokamak geometry with ε = 0.1 in the absence of the Shafranov shift, plasma elongation and triangularity (equilibrium density and temperature gradients are L n /L Tj = 1 with Φ ′ eqm = 0). The threshold we derive is the result of the radial shift of drift islands described by the S function, equation (21), and, in particular, the pressure gradient restoration across the magnetic island O-point at w ∼ ρ ϑi . As discussed in section 5, the latter mainly arises from the behaviour of the σ-dependent part of the ion distribution function, σ g (0,0),σ i , at small w. Therefore, we emphasise that this threshold physics is related to passing particle dynamics, and not the finite banana width effects of the trapped particles.

Summary and conclusions
To summarise, a new drift kinetic theory of magnetic islands, valid for w ∼ ρ ϑi , has been developed in a low collisionality plasma. The electron/ion distribution function is found to be flattened across the so called drift islands, which are radially shifted by a value, proportional to ρ ϑe/i /w, compared to the magnetic island. This, in turn, results in a density gradient being sustained across a magnetic island of width comparable to the ion poloidal Larmor radius, w ≲ ρ ϑi . For ions, their finite 10 As shown in appendix A, the (1/B)∂B/∂ψ term is ε times smaller than the rest of the terms kept in equations (14) and (15). Therefore, L −1 B = 0 provides results, close to the case when the actual ϑ dependence is maintained in LB, which we can see from figure 8. density gradient at the centre of the magnetic island is the consequence of the effect of the drift islands, i.e. the relatively large radial shift in equation (21). While the electron radial shift is small since ρ ϑe ≪ ρ ϑi , their density gradient is still present across the magnetic island O-point due to the effect of the electrostatic potential that is generated to ensure plasma quasi-neutrality.
This gradient is found to suppress the drive for NTMs when w is small, providing the critical magnetic island width, below which NTMs decay. We note that this physics is associated with the passing particle dynamics and not the effects of the finite banana orbit width. In this sense, the relevant parameter for the NTM threshold would be the ion poloidal Larmor radius. However, we find that in addition to the linear dependence of w c on ρ ϑi , w c also scales with the inverse aspect ratio as ε 1/2 . This indicates that the threshold is also influenced by the trapped particles-this is physics that we do not yet fully understand. For the small inverse aspect ratio circular poloidal cross section tokamak, it scales as w c = 0.45ρ ϑi at ε = 0.1 in the island rest frame. This provides w c = 1.41ρ bi , if written in terms of the ion banana orbit width. This result is in quantitative agreement with the previous experimental observations, e.g. [33].
If we consider contours of the total density, obtained with the self-consistent electrostatic potential and differentiated with respect to ψ, we notice that these contours in the (ψ, ξ) plane would reproduce the contours of the flow parallel to the field lines. This changes abruptly with ρ ϑi /w, as can be seen from figures 4(a)-(d). Therefore, comparing the total density profiles that contain the magnetic island perturbation against the experimental density profiles for a large magnetic island and for a small island near the marginal stability might be one of the steps to check the validity of the drift island concept as an explanation of the threshold, w c .
The theory presented here provides a prediction for the magnetic island threshold obtained in the conventional tokamak limit in the island rest frame. However, as mentioned above, the theory presented in this paper is subject to certain limitations. One of them is the boundary layer around the drift island separatrix. In this layer, the transport terms will become comparable to the parallel streaming, i.e. the left hand side of equation (20), ν ii/e ∂ 2 /∂S 2 ∼ A ∂/∂ξ| S . This layer can be treated in a similar way to the collisional dissipative layer at the trapped-passing boundary, which we considered in section 5 by exploiting its thinness. However, as the drift islands are shifted radially relative to the magnetic island by an amount that varies with V ∥ and is proportional to ρ ϑj , this might result in a broader layer around the magnetic island for the electrostatic potential, being spread over a width comparable to ρ ϑj . It is this layer in the potential that influences the polarisation current when ω is not equal to zero. We will present this in future publications, extending existing theories of the separatrix layer [21,22,29,30]. A special challenge is where both boundary layers overlap, A ∂/∂ξ| S ∼ ν ii/e ∂ 2 /∂λ 2 ∼ ν ii/e ∂ 2 /∂S 2 . However, since this area is very narrow in the low collisionality plasma, we can fix all the coefficients in the full equation, equation (14), and thus find the solution even in the presence of the electrostatic potential in a way similar to section 5. The obvious advantage of this method, employed in RDK-NTM, is that it allows us to obtain the solution with sufficient resolution in these narrow boundary regions. Another option to deal with the separatrix layer in S space is to solve equation (14) for the full 4D distribution function as a function of p φ , ξ and λ at each V. This is the DK-NTM approach, presented in appendix C. However, as mentioned above, it is then difficult to efficiently resolve both layers simultaneously and there are numerical limits to how small we can take ν j or ρ ϑi /w. For the collision operator employed in this paper, which captures aspects of neoclassical transport, the thinness of the boundary layer that surrounds the drift island is provided by low ν ii/e . More generally, the cross-field transport will be driven by turbulence, and then the layer physics will be much more complicated to capture, especially as the island and the layer length scales become comparable to turbulent eddy size-this will become a challenging multiscale numerical problem. Both the 3D and 4D approaches have their benefits and drawbacks, and both require extra development to refine the treatment of the separatrix layers. We highlight that while the proper treatment of the separatrix boundary layer is necessary for calculating accurately the layer contribution to the polarisation current, this does not influence the bootstrap current. Thus, the result presented here for the NTM threshold when the island is at rest in the plasma E × B reference frame is expected to be robust.
The second problem that is beyond the scope of this paper is the calculation of the island rotation frequency [39]. This is associated with the dissipation processes in the plasma. In [39] it is found from the solution in the dissipation layer around the trapped-passing boundary, provided it is the main source of dissipation. Potentially, there might be an extra contribution from the separatrix layer and thus is to be further investigated. While the bootstrap current perturbation described here is not expected to depend on the island propagation frequency, the polarisation current will, and this could influence the threshold.
Another problem to be addressed in the future is the curvature contribution to the island evolution. The drift kinetic equation presented in this paper is correct to O(ε 2 ∆g j ), which is already sufficient to find ∆ cur . The more accurate calculation, though, would also require the O(ε 2 ∆ 2 g j ) terms. We highlight that while the effects of the plasma shaping are weak in conventional devices, they would be more important for spherical tokamaks and potentially further enhance the stabilising mechanisms that provide a threshold.

Data accessibility statement
The DK-and RDK-NTM numerical data used for benchmarking is available at https://doi.org/10.5281/zenodo.4294044.

Acknowledgments
This work was supported by the UK Engineering and Physical Sciences Research Council, Grant No. EP/N009363/1. Numerical calculations were performed using the ARCHER computing service through the Plasma HEC Consortium EPSRC Grant No. EP/R029148/1. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

∂R ∂ψ
with I ′ ∼ R 2 p ′ /I and hence I ′ /R ∼ β/r 2 . Thus, in a low beta limit, we find, with R = R 0 (1 + εcosϑ) taken for the model major radius.
Since ∂/∂λ| ψ and ⟨...⟩ pφ ϑ are not commutative, we transform the pitch angle derivatives, using . Equation (B3) allows us to write: Then applying equation (B4) to switch fromp φ to S, we obtain: Here The first two terms in equation (B6) will provide the pitch angle scattering in S space, while the other three terms are responsible for transport across surfaces of fixed S. Substituting equation (B6) into the above collision operator providesĈ j and hence the collisional constraint in S space, equation (23). After multiplying both sides of equation (23)     Same as Figure C2, except that λ is close to the trapped-passing boundary, λc = 0.91. by the DK-NTM model; this is assumed to be numerical, but will be explored further in future. Since it is the density/temperature gradient that is responsible for the perturbation in the bootstrap current, in figure C4 we show the density moments, differentiated with respect to p φ and plotted against p φ to magnify any differences. The largest discrepancies are at the largest ρ ϑi /w. We note from these comparisons that the (R)DK-NTM distribution functions agree reasonably well even around the S island separatrix, where the ∂ 2 /∂S 2 term from the collision operator competes with the parallel streaming.
In figures C2-C4 we confirm that the impact of the drift islands on the density gradient inside the magnetic island increases with ρ ϑi as the radial shift of drift islands relative to the magnetic island becomes larger. At the same time, the effect of the radial shift seems to become less crucial as λ approaches the collisional dissipation layer, e.g. the ion density moment has a partially flat region across the island at λ = 0.83 even at a relatively large ρ ϑi /w = 0.35 or 1. This might be explained by the fact that the passing solution close to λ = λ c is influenced by the trapped branch via the matching in the dissipation layer and vice versa: drift islands do not exist in the trapped region, but the magnetic island is still present and provides the flat solution in the vicinity of the O-point. Nevertheless, since the fraction of trapped particles is small, their contribution to the total density is insufficient to enforce the conventional flat profile across the magnetic island at large ρ ϑi w −1 .
In figure C5 we show the flow moments, (1/2) σ σg (0),σ i , in ψ space. As mentioned above, since the trapped particle solution is independent of σ at fixed p φ , their contribution to the parallel flow is small, as we can see from figure C5(a). Both the DK-and RDK-NTM solutions demonstrate agreement with the drift island concept described above in the vicinity of the magnetic island and match the equilibrium neoclassical flows far from the island. At small ρ ϑi , e.g. ρ ϑi /w = 0.05, the DK-NTM solution has oscillations in the flow as ξ approaches the X-point as shown in figures C6(c) and C7(b). This is potentially numerical, noting that the diffusion term in p φ or S space scales with ρ ϑi and will become difficult for the 4D DK-NTM model to resolve at small ρ ϑi /w. This is true even at a relatively large ν * i = 10 −2 , e.g. see figures C2(g) and C3(g). This also explains why the best matching results are obtained at ρ ϑi ≲ w, while achieving smaller ρ ϑi /w becomes increasingly challenging for the 4D DK-NTM code as the transport terms become smaller. A similar discrepancy is found in collisionality. Indeed, the pitch angle scattering outside the dissipative layer is small and is dominated by the free streaming. A lower ν * j will decrease the collision operator further, making it difficult for the DK-NTM to resolve the collisional terms. In contrast, the RDK-NTM model requires low collisionality for particles to follow the S streamlines and to implement matching through the collisional boundary layer. We find, ν * i ∼ (10 −3 , 10 −4 ) provides a narrow window where the validity regions of both solutions overlap.
To investigate the basis for the threshold, in figures C7(a), (b) and C8(a), (b) we integrate the flow moments over velocity space to provide the ion flow for both the DK-and RDK-NTM models for small and large ρ ϑi . At ρ ϑi /w = 0.05 (figures C7(a) and (b)), both the DK-and RDK-NTM reproduce the results of the conventional model for large magnetic islands with a zero flow profile across the island O-point and matching the neoclassical equilibrium far from the magnetic island. As discussed above, the flat region within the magnetic island decreases with increasing ρ ϑi /w, as can be seen by comparing Figures 4, C8(a) and (b). At the same time, the peak region of the parallel flow right outside the separatrix becomes wider and is compensated by the contribution to the parallel flow that originated inside the island. As we can see from figures C8(a) and (b), this balancing contribution is mainly localised within the island region according to the DK-NTM predictions (blue area inside the island in figure C8(b)), while in the RDK-NTM solution it is concentrated around the X-point (blue areas closer to the X-point in figure C8(a)).
The origin of the negative DK-NTM flow in figure C8(b) inside the magnetic island is shown in figure C10. In contrast to the RDK-NTM solution, where σ = ±1 drift islands are located at the same level with the magnetic island, in accordance with equation (21), the DK-NTM distribution functions (and hence drift islands) for σ = ±1, along with the radial shift, are also vertically shifted relative to each other and to the magnetic island. The vertical shift of the distribution function is denoted by g 0 in figure C10. As we see from figures C10(b)-(d), DK-NTM predicts that g 0 is not constant but is a function of λ at fixed p φ and/or ψ. At fixed p φ and λ, the vertical shift is a function of ξ as well. Its origin in DK-NTM is to be investigated further as part of the future work. Note, as the DK-NTM solutions are vertically shifted by the same value, this vertical shift does not influence the density moments and thus they agree with the RDK-NTM solutions. Now let us compare the (R)DK-NTM solution against the analytic solution found in [18]. The curvature of the distribution function in the vicinity of the island separatrix is determined by the diffusion terms that arise from transforming from ψ to S/p φ in the pitch angle scattering collision operator. These transport terms are proportional to ∂ k /∂S k λ,ξ or ∂ k /∂p k φ λ,ξ (k = 1, 2), respectively. In figure C9(c) we plot the leading order ion distribution function, differentiated with respect to y, against y, obtained with the full collision operator. For RDK-NTM this is the collision operator presented in appendix B.
In figure C9(d) we plot the analytic result of [18] valid in the limit of large w (ρ ϑi /w = 0.05 is chosen). This is to be denoted by H96. H96 is derived from a model diffusion of the form Γ ψ = −D∂n/∂ψ, where Γ ψ is the particle flux in the radial direction and D is the diffusion coefficient that has been assumed to be a slowly varying function across the magnetic island Opoint. In figure C9(d) we also show the test RDK-NTM solution with the model radial diffusion weighted by √ S + cos ξ (green markers) 12 . The latter reproduces well the H96 solution at large w outside the island. Figures C9(c) and (d) show that keeping the actual S diffusion is important; otherwise, a      significant fraction of the perturbation will be removed right outside the separatrix. A consequence of this is that our quantitative results may be sensitive to turbulent diffusion-a possibility to explore in the future.
A potential source of discrepancy between the DK-and RDK-NTM models is associated with a narrow boundary layer around the drift island separatrix. While the DK-and RDK-NTM g (0) i seem to agree well in the vicinity of the separatrix, figure C9(c). shows that there is some difference in their y derivatives. In the vicinity of the drift island separatrix, there is a region where S derivatives drive a large diffusion which can be comparable to parallel streaming, ν ii/e ∂ 2 /∂S 2 ∼ A ∂/∂ξ| S . This would invalidate the perturbative treatment of collisions in the RDK-NTM solver. This region is then to be treated in a way similar to the dissipative layer solution in λ space discussed above. This separatrix layer has been addressed to some extent in [21,22,29,30], but without a complete treatment in toroidal geometry. Future extensions of the RDK-NTM model will address this, but it is not expected to influence the calculation of the bootstrap current and thus the result for the RDK-NTM magnetic island threshold presented in section 6. However, the difference shown in figure C9(c) might be crucial in calculating the polarisation current, which at this stage remains part of our future work (noting that here we consider stationary islands in the plasma E × B rest frame, for which the polarisation current is not expected to be important).