Finite ion orbit width effect on the neoclassical tearing mode threshold in a tokamak plasma

A new drift kinetic theory for the response of ions to small magnetic islands in toroidal plasma is presented. Islands whose width w is comparable to the ion poloidal Larmor radius are considered, expanding the ion response solution in terms of , where r is the minor radius. In this limit, the ion distribution can be represented as a function of toroidal canonical momentum, . With effects of grad-B and curvature drifts taken into account, the ion distribution function is a constant on a ‘drift island’ structure, which is identical to the magnetic island but radially shifted by . The distribution is then flattened across the drift island, rather than the magnetic island. For small islands , the pressure gradient is maintained across the magnetic island, suppressing the bootstrap current drive for the neoclassical tearing mode (NTM) growth. As , the ions are largely unperturbed. However, the electrons respond to the electrostatic potential required for quasi-neutrality and this provides a stabilizing contribution to the NTM evolution. This gives a new physical understanding of the NTM threshold mechanism, with implications for the design of NTM control systems for future tokamaks such as ITER.


Introduction
A neoclassical tearing mode (NTM) is an example of a resis tive magnetohydrodynamic (MHD) instability present in a tokamak plasma. NTMs are characterized by the evolution of a magnetic island chain, which arises from a filamentation of the component of the plasma current density parallel to the magnetic field lines on a rational surface. The change in magn etic topology enhances the transport of particles and heat across the island width, as they stream along the perturbed magnetic field lines. Consequently, in the absence of particle/ heat sources and particle drifts, the radial pressure gradient is flattened across the island width. This not only degrades the plasma confinement by lowering the pressure in the core of the plasma, but a large NTM can trigger a disruption. Major disruptions can do significant damage in future larger toka maks such as ITER; it is therefore crucial to control or avoid NTMs.
The evolution of a magnetic island can be described by the modified Rutherford theory. According to the original theory [1], the stability of an island depends on the tearing param eter, ∆ , which is a measure of the free energy available in the current density for magnetic reconnection. In the modified theory incorporating toroidal geometry, a localized perturba tion in the parallel current profile is also considered. One such contribution comes from the bootstrap current, which tends

Nuclear Fusion
Finite ion orbit width effect on the neoclassical tearing mode threshold in a tokamak plasma K. Imada 1 , H.R. Wilson 1,2 , J.W. Connor 2 , A.V. Dudkovskaia 1 and P. Hill 1 to be removed from the island region where the pressure gra dient is flattened. This perturbation in the bootstrap current typically reinforces the original filamentation in the current density profile, which enhances the island growth [2,3]. This enhanced island growth mechanism, provided by the per turbed bootstrap current, is what characterizes a neoclassical tearing mode. Because this bootstrap current contribution (which we label ∆ bs ) scales as 1/w, where w is the island half width, this theory predicts that all seed islands, however small, would grow to a large saturated island (assuming a negative ∆ ). If this were the only contribution to the island evolution, then it would pose significant challenges for achieving fusion in a tokamak, as all plasma discharges would suffer from such tearing instabilities. However, experimental observations point to the existence of a threshold effect [4,5], whereby a sufficiently small seed island heals itself and shrinks away. One of the possible origins of this threshold is the finite radial transport effect [6][7][8]. This effect partially restores the pres sure gradient that is flattened across the island width, which reduces the bootstrap current drive for the island growth. Then, the growth of an island with a width comparable to or smaller than w χ /r = ( nL q /r) −1/2 (χ ⊥ /χ ) 1/4 is suppressed, and can be stabilised for negative ∆ . Here, = r/R is the inverse aspect ratio, r and R are minor and major radii of the torus respectively, n is the toroidal mode number, L q is the magnetic shear length scale and χ ⊥, are perpendicular and parallel thermal conductivities, respectively.
Another candidate for producing a threshold is the polari zation current, which is induced when the island chain prop agates through the plasma with a characteristic frequency, ω. In slab geometry, the finite Larmor radius (FLR) effect sets the length scale where the polarisation current becomes relevant. When w is comparable to the ion Larmor radius, ρ Li , the dif ference in the responses of ions and electrons to the rotating island gives rise to an electrostatic potential. The difference in the gyroaveraged E × B drifts then leads to the polarization current, which in turn generates a parallel return current, δJ , to ensure ∇ · J = 0. This δJ then contributes to the magnetic island evolution through the modified Rutherford equation. Previous works [9][10][11][12] have shown that there exists a narrow boundary layer current in the vicinity of the island separatrix, whose contribution to the island evolution is opposite but com parable in magnitude to that from outside the boundary layer. Whether or not the polarization current contribution is stabi lizing depends on the relative sizes of the two contributions.
In toroidal geometry, the combination of gradB and curva ture drifts causes orbits of passing particles (those completing full poloidal orbits) to stray from a reference flux surface by a distance of O( ρ θ ), where ρ θ is the poloidal Larmor radius. In addition, because of the variation in the magnetic field strength, a fraction of particles are trapped on the outboard side of the tokamak. They execute closed bananashaped orbits, whose width is given by ρ b ∼ √ ρ θ . The trapped ions experience a different orbitaveraged E × B drift to the trapped electrons, which results in a net current: the neoclas sical polarization current [13][14][15][16]. Its contribution to the island evolution (which we label ∆ pol ) can be substantial when w is comparable to the ion banana width, ρ bi . Previous works [14,17] based on drift kinetic theory in toroidal geometry have shown that ∆ pol ∝ 1/w 3 , when ρ bi w. If this contrib ution is stabilizing, then it could heal small seed islands, thus providing the threshold.
One way of controlling NTMs is to use electroncyclotron current drive (ECCD) [18,19]. An ECCD system can suppress the NTM growth by driving a plasma current in the vicinity of the island, effectively replacing the missing bootstrap cur rent in the region. If the island width can be reduced to below the threshold width, then the island will shrink away and the NTM will be successfully suppressed. The effectiveness of the ECCD has been demonstrated in a number of tokamaks [20][21][22][23], and it is currently the favoured method of control ling NTMs in ITER [24]. However, the power consumption of the ECCD system is rather high. It is therefore crucial to use it as efficiently as possible, if we are to achieve a high fusion Qfactor in future experiments [25]. Here, Q is the ratio of fusion power to the heating power injected into the plasma. This is why the understanding of the NTM threshold physics, including the predictive capability for the threshold island width, is so crucial.
Taking into account these contributions, the modified Rutherford equation describing the island evolution can be represented in the following form: where τ R is the resistive time scale, β θ = 2µ 0 p/B 2 θ , p is the plasma pressure, B θ is the poloidal component of the magnetic field, L q = q(dq/dr) −1 , q is the safety factor, L p = −p(dp/dr) −1 , g( , ν ii ) describes the collision frequency dependence [17], and a 1,2 are numerical constants. The terms in a 1 and a 2 correspond to the bootstrap current (∆ bs ) and polarization current (∆ pol ) contributions, respectively, while ∆ ECCD describes the impact the ECCD system has on the island evolution. As mentioned above, an efficient deployment of the control system requires a better understanding of the threshold physics. In particular, the ability to quantitatively predict the threshold island width w c (for which dw/dt = 0), below which the island shrinks away, is essential. Existing theory relies on the assumption ρ bi w, i.e. valid in the limit of large island widths compared to the ion banana width. However, observations [26] show that the threshold island width is often comparable in size to the ion banana width; precisely the regime where this assumption breaks down. A new theory is therefore required to accurately determine the relative sizes of the ∆ bs , ∆ pol and ∆ contrib utions, including their dependence on the curvature and finite particle orbit width effects. This will allow us to quantitatively predict the threshold width for ITER. This paper focuses on the effect of finite ion orbit width on the bootstrap current contribution to the magnetic island evo lution, by considering a stationary island relative to the plasma. We extend the existing drift kinetic theory [14] to describe the ion response to the island perturbation. Crucially we relax the small bananawidth assumption, and consider magnetic islands whose widths are comparable to the ion banana width: w ∼ ρ bi . We consider a small magnetic island w r (valid for an island with a width close the threshold width), which allows us to treat the plasma as toroidally symmetric to leading order. Then, because of the finite orbit width effect the ion distribu tion function is no longer a function of poloidal magnetic flux ψ, but of toroidal canonical angular momentum: which is conserved along the orbits in an axisymmetric plasma. Here, ψ s is the poloidal flux at the rational surface where the island is located, I = RB φ , B φ is the toroidal component of the magnetic field, v is the component of the particle velocity along the field lines and ω ci is the ion gyrofrequency. For electrons, the assumption ρ be w is still valid, which allows us to use the existing analytic solutions for the electron distribution function of [14]. However, in the regime ρ be w ∼ ρ bi , we anticipate a notable difference in the electron and ion distribution functions, if we neglect the electrostatic potential, Φ. It is therefore important to cal culate Φ selfconsistently from quasineutrality, which we incorporate into our analysis. From the particle responses, we determine the full contribution of the localized current perturbation (which includes the perturbed bootstrap current) to the island evolution, including those from inside the island and the separatrix layer (previously assumed to be zero, with perfectly flat pressure gradient inside the island separatrix). Earlier works considered the limit w ∼ ρ bi by employing a particleincell (PIC) simulation to solve the drift kinetic equation [27], or approached the problem analytically by focusing on the contribution from the passing particles only [28]. They found that an ion density gradient is supported across the island, but did not address the consequences this has for quasineutrality and the electron response. Our ana lytic approach reveals the physics explanation for the density gradient, and provides a new threshold physics effect that results from the electron response. This paper is organized as follows. In section 2, we intro duce the perturbed magnetic geometry and drift kinetic equa tion, followed by section 3 outlining the analytic electron response, as derived in [14]. Electron flow depends on the ion counterpart through the model collision operator, and hence is worthwhile revisiting here. In section 4, we derive the orbitaveraged equation describing the ion response and first consider the solution in the collisionless limit, where we introduce the concept of 'drift islands' in shifted flux space (section 4.1). This aids the understanding of solutions of the full equations, including the collisional effects. This is followed by section 4.2, describing the analysis of the full solution for the ion distribution function, perturbed density and parallel flow profile. Finally, in section 5 we determine the layer current contribution to the island evolution, ∆ loc , for a range of values of ρ θi and w. Conclusions are drawn in section 6.

Magnetic island geometry and drift kinetic equation
We consider a large aspect ratio, circular crosssection tokamak, neglecting the Shafranov shift. Then, in the orthog onal coordinate system ∇φ × ∇ψ = rB θ ∇θ, where ψ is the poloidal magnetic flux, θ is the poloidal angle and φ is the toroidal angle, the equilibrium magnetic field is given by: Here, I(ψ) = RB φ and B θ and B φ are the poloidal and toroidal components of the magnetic field respectively. A magnetic perturbation of the form satisfying Maxwell's equation is introduced: where b 0 = B 0 /B 0 is the unit vector in the direction of the equilibrium field lines, and the parallel vector potential takes the form: assuming a single dominant helicity perturbation. Here, ξ is the helical angle in the island rest frame: where m is the poloidal mode number, q s = m/n is the safety factor at the rational surface where the island is located (all quantities with subscript s are those evaluated at the rational surface, unless otherwise indicated), and n is the toroidal mode number. ψ = (w 2 ψ /4)(q s /q s ) describes the perturbation amplitude (the prime denoting a differential with respect to ψ), and w ψ is the island halfwidth in ψspace. It is also con venient to introduce a perturbed flux function Ω that describes the magnetic island geometry, which satisfies: B · ∇Ω = 0. That is, the perturbed magnetic field lines lie in surfaces of constant Ω. The form of Ω is then given by: with Ω = 1 defining the island separatrix. Working in the island rest frame, we consider a steady state ion response to the magnetic island perturbation, using the drift kinetic equation: for a particle species j . Here, denotes a comp onent parallel to the magnetic field lines, is the combination of gradB and curvature drifts, ω cj = e j B/m j and e j and m j are the par ticle charge and mass respectively. Φ is the perturbed electro static potential to be determined from quasineutrality, and C j is the momentumconserving model collision operator [29]. Likelike particle and electron-ion collision operators are respectively given by: where the λ differentials are taken at fixed ψ, is the deflection In equation (9), we have also introduced: which is required for momentum conservation, and In equation (8), spatial derivatives are taken at constant kinetic energy, E = v 2 /2, and magnetic moment, µ = v 2 ⊥ /2B, where ⊥ denotes a component perpendicular to magnetic field lines. Working in (v, λ) velocity coordinates, where λ = µ/E is the pitch angle, the velocity space integral is: where σ is the sign of parallel velocity, v = σv √ 1 − λB.

Electron response
In order to calculate the perturbed bootstrap current and its contribution to the magnetic island, we require both the ion and electron parallel flows. For electrons, we can exploit the fact that the electron poloidal Larmor radius is usually small compared to the magnetic island width, even when the ion poloidal Larmor radius is comparable to the island width (i.e. ρ θe ρ θi ∼ w). Then, we are justified in using the analytic results for the electron response, derived in [14]. As can be seen from equation (10), the electron parallel flow is depen dent on the ion counterpart through momentum conservation, so the two species are coupled.
We seek a Maxwellian solution for the electrons, expanding the distribution function about the rational surface in the limit of a small island w r: where −e is the electron charge, g e describes the perturbation in the electron distribution function, is the Maxwellian for species j and n 0 is the equilibrium den sity. Introducing two small parameters δ e = ρ θe /w and ∆/r, we expand the perturbation term: and determine the leadingorder response of electrons to the island. Then, considering relevant order contributions to the drift kinetic equation (8), the solution for the electron distribu where h(Ω) describes the perturbed radial density profile in the vicinity of the magnetic island. In the absence of the drift effects, the gradient would be completely flat inside the island separatrix, Ω < 1 (note the Heaviside function in equation (18)). h e in equation (17) is determined from a constraint equa tion derived from the O(δ e ∆) contribution to the drift kinetic equation, which takes the form: where g (1,0) e is the last two terms of equation (17): Equation (19) can be solved analytically in the collisionless limit, when the term in the collision operator is negligible. Solving (19) for h e yields: where · · · θ denotes an integral over a period in θ for passing particles and between the bounce points, multiplying the result by σ and summing over σ for trapped particles. · · · Ω denotes a flux surfaceaverage: H e (Ω) is a free function satisfying: Solving for H e using the result for h e from equation (20), and integrating over velocity space, the result for the flux surface averaged electron parallel flow is: where f t and f c are trapped and passing particle fractions respectively and k −1.173 [30]. (Hydrogenic, quasineutral plasma is assumed for simplicity.) We interpret the flux sur faceaveraged parallel flow as the component of the flow that gives rise to the bootstrap current (see equation (35) later). This result for the electron parallel flow will be used in sec tion 5 to calculate the perturbed bootstrap current, as well as its contribution to the island evolution.

Ion response
For ions, we relax the assumption of small ion poloidal Larmor radius relative to the island width and consider an arbitrary ratio ρ θi /w, which is the critical difference from the treatment of electrons in the previous section and past analytic works [14,17]. Taylorexpanding the Maxwellian about the rational surface where the island is located (ψ = ψ s ) in the small island limit w r, we seek a solution to the ion distri bution function in the vicinity of the magnetic island: Using the parameter ∆ = w/r 1, this time we expand the perturbation in the ion distribution function in terms of ∆ only, retaining the ordering ρ θi ∼ w: When ρ θi ∼ w, both parallel streaming and magnetic drift dominate the ion response. Then, the leading order contrib utions to the drift kinetic equation (8) where ω * i = mcT i n /Zeqn is the ion diamagnetic frequency, In the limit of a small island perturbation, the toroidal symmetry is approximately conserved to the leading order in ∆. Then the toroidal canonical momentum (2) is a conserved quantity along particle orbits, which we can utilize as a radial coordinate in place of ψ. As we shall see, this allows us to eliminate one of the spatial coordinates, θ. Thus, transforming from radial variable ψ to p φ , equation (25) simplifies to: v Rq which can straightforwardly be integrated to yield: where the free function, h 0 , is determined from the next order equation. This allows us to write the total distribution function as a function of p φ instead of ψ: whereḠ The physical meaning of this is that the distribution function is a constant on the orbits the particles freestream along, rather than being a flux surface quantity. These orbits are described by p φ = constant (i.e. the standard neoclassical orbits for a toroidally symmetric system). h 0 then describes the modifi cation to the equilibrium profile of the distribution function, when the island perturbation is introduced. We see later that the finite orbit width effect has a profound impact on the per turbed density profile in the vicinity of a magnetic island. We now proceed to the O(∆) contribution to the drift kinetic equation (8): v Rq where C ii is given by equation (9). To eliminate the term in g 1 , we take the average of equation (30) along particle orbits, at fixed p φ . For passing particles, this is achieved by multi plying equation (30) by Rq/v and integrating over a period in θ at fixed p φ , making use of the periodicity in g 1 . For trapped particles, the distribution functions at the bounce points satisfy: g 1 (σ = +1, θ b = ±1) = g 1 (σ = −1, θ b = ±1), by conservation of particles. Thus, the term in g 1 for trapped particles can be eliminated by multiplying equation (30) by Rq/|v |, summing over σ and then integrating with respect to θ between the bounce points. The result is a particle orbit averaged equation for h 0 (through Ḡ 0 ): where ψ = p φ + Iv (θ)/ω c . Note that these θ integrations on the ion equation are performed at fixed p φ , while those for electrons were at fixed ψ.
Since we are concerned with a Maxwellian solution in the vicinity of the island (the equilibrium profile is assumed far away from the island), we consider, for simplicity, a large aspect tokamak with B = 1 − cos θ. The first line of equa tion (31) can then be expanded about the rational surface. Then, introducing normalized quantities: we obtain the dimensionless equation Here, Ĉ ii = (Rq/v thi )C ii and y c = 1 corresponds to the trapped/passing boundary in the pitch angle space. Before solving equation (32) in full, in the next section we consider the form of the solution in the collisionless limit.

Collisionless limit
In section 2 we introduced the perturbed flux function Ω describing the magnetic island geometry (see equation (7)). The perturbed magnetic field lines lie in the surfaces of con stant Ω. In our present analysis, we introduce a new set of surfaces defined by S: Note that, for Φ = 0 and y < y c (passing particles) the con stant S surfaces are identical to the constant Ω surfaces, but shifted radially by an amount proportional to ρ θi . This shift can be attributed to the term in ρ θiωD in equation (33) and the second term of p φ (see equation (2)). Working with S as the new 'radial' coordinate, we can further simplify equation (32) for Ḡ 0 , which now takes the form: where it should be noted that the differential with respect to ξ is now taken at fixed S. This illustrates that the streamlines lie in surfaces of constant S, not constant Ω. They differ because of the particle orbits-the radial shift of the 'drift island' of the constant S contours relative to the magnetic island is due to the gradB and curvature drifts, while the term in Φ arises from E × B drifts. A similar radially shifted structure has been found for the ion flux, which originates from the addi tional guiding centre drift caused by the electrostatic potential perturbation in the vicinity of the magnetic island [31].
In figure 1  In the absence of the electrostatic potential term, this shift is equal and opposite for the v < 0 case. We call this shifted island structure in S the 'drift island', whose physical conse quence is paramount. In the low collision frequency limit, where we may assume that the term on the right hand side of equation (34) becomes O(ν ii Rq/v thi ), 1 smaller, it can be shown straightforwardly that the solution for Ḡ 0 is a function S: , v). That is, the constant S surfaces are also constant Ḡ 0 surfaces, meaning the region over which the radial gradient of the distribution function is flattened is also shifted to coincide with the drift islands rather than the magnetic island. Because the shifts in the distribution function are in opposite directions for σ = +1 and σ = −1, when the perturbed density is con structed by a velocity integral summing over σ, a substantial radial gradient is supported inside the island when w ∼ ρ θi . For a large island, w ρ θi , the shifts are relatively small and the density moment is then approximately flattened across the island, as expected. This is what we call the finite orbit width effect, which is distinct from the well known radial transport effect [6]. It lies at the heart of the restored density gradient found in the PIC simulations of [27]. In the next section, we demonstrate this effect with the full solution to equation (32), as well as discuss the consequences for the parallel flow profile and hence the impact on the current and island evolution.

Full solution for ion response
In figure 2, we present the colour contour plot of the full per turbed ion distribution function f i in the x − ξ plane obtained by solving equation (34) numerically for ν * = 0.01. The Boltzmann factor with the perturbed electrostatic potential is included in the plot. This has been determined via quasi neutrality, using the electron response derived in the previous section. The solid lines are contours of constant S, while the dashed line indicates the location of the magnetic island sepa ratrix. The plot clearly shows that the colour contours of f i are wellaligned with the contour lines of S, and they are radially shifted relative to the island separatrix. This indicates that f i is indeed a function of S to leading order, and the radial shift of the profile is O(ρ θi ), as expected. As described in the previous subsection, the equal and opposite shifts for σ = ±1 have a significant consequence for the radial density profile. As shown in figure 3, flattening of the density gradient inside the magnetic island is wellpreserved for ρ θi w but is almost absent for ρ θi ∼ w. The restoration of the density gradient across the magnetic island is precisely the result of the shift in the drift islands; because the flat spots in the shifted distribu tion functions for σ = ±1 no longer align when w ∼ ρ θi , the summation over σ causes the gradient to be maintained across the island. Specifically, the σ = +1 solution for f i has a gra dient where the σ = −1 solution is flattened, and vice versa (see figure 4). On the other hand, if ρ θi w, then the flat regions for σ = ±1 do align to a large extent and the density gradient is flattened inside the magnetic island, as expected.    For electrons, the strong parallel flow tends to keep the density flattened across the magnetic island width, even for small islands (i.e. ρ θi ∼ w, but ρ θe w). However, the elec tron distribution function depends on the electrostatic poten tial as well, in such a way as to satisfy quasineutrality. This selfconsistent potential influences both the electrons and ions. Therefore the full density, including the Boltzmann factor, takes the form given in figure 3 for both ions and electrons. This physics has consequences for the structure of the electro static potential (figures 5 and 6). When ρ θi w, the potential is constant on the perturbed flux surfaces, as expected from previous theories. However, when ρ θi ∼ w, this is no longer the case. Furthermore, the region inside the island retains a substantial potential gradient, consistent with the picture described above. The same is true for the ion parallel flow profile, as shown in figure 7. For large islands, the flow is a perturbed flux quantity, with a well defined boundary layer flow in the vicinity of the island separatrix. Conversely, for a small island the flow is no longer constant on the flux surfaces, and the boundary layer structure is completely lost.

Contributions to island evolution
We now have all the elements to consider the contribution to the island evolution originating from the perturbed current, localized in the vicinity of the rational surface, ∆ loc . We dis tinguish the bootstrap current contribution from other sources such as the neoclassical polarization current by defining it as the component of the total current that is constant on the per turbed flux surfaces. The bootstrap current can then be pro jected out by: where ... Ω is given by equation (21). In previous analytic works [14,17], this J Ω is interpreted as solely consisting of the flux surface average of the perturbed bootstrap current.
Here, however, when the island is small (i.e. w ρ θi ) there is an additional contribution that cannot be explained purely in terms of the standard bootstrap current picture, as described later. Therefore we refer to it as the localized current perturba tion, with its contribution to the island evolution labelled as ∆ loc . It is calculated from the dispersion relation derived from Ampère's law: Figure 8 shows the results for ∆ loc normalised to β θ = 2µ 0 p/B 2 θ as a function of w for a range of values of ρ θi . For large w ρ θi , ∆ loc tends to the asymptotic value (lim ρ θi /w → 0) expected from previous analytic theories for the bootstrap drive [2,14], which is represented by the dashed line. However, for small island widths approaching the size of ρ θi , we see that the impact of the shifted drift islands is to reduce the bootstrap drive. For even smaller island widths, ∆ loc becomes negative. This is a rather remarkable result, as  it means that the effect of the current perturbation is to heal the island and therefore represents new threshold physics that cannot be explained by a reduced bootstrap drive alone. For larger ρ θi , the peak value in ∆ loc decreases substantially, hence suppressing the bootstrap drive for the island growth. The critical island width, w c , where ∆ loc passes through zero, increases linearly with ρ θi : it can be fitted by w c 2.76ρ θi (see figure 9). Experimental observations support this linear relationship [26], though the coefficient we derive is some what larger than the result obtained from experiments.
We now consider the physics underpinning the stabiliza tion of small islands, w ρ θi . Figure 10 shows the plots of ρ θi × ∆ loc /β θ versus w/ρ θi , with separate ion and electron contributions. It is clear that all cases with different ρ θi /r values condense onto a universal set of curves for both the ion and elec tron contributions. This is a consequence of the parallel flows being proportional to ρ θi,e , as predicted by analytic neoclassical theory. An important point to address is that, as w → 0, the ion contribution to ∆ loc tends to zero. This is consistent with the density gradient (and therefore bootstrap current) being unper turbed in this limit (as found in the PIC simulations of [27]). Indeed, we expect that when the island width is much less than the ion banana width, the ions will average over the perturbed electromagnetic fields associated with the island. Electrons still respond to the perturbed fields, and we see from figure 10 that   [14] for the bootstrap current contribution, for which ∆ bs ∝ 1/w.  it is their response that provides the stabilizing contribution. We postulate that it is the response of the electrons to the electro static potential, required for quasineutrality, that creates the stabilizing contribution to the current density.

Conclusion
We have presented a new drift kinetic theory for the response of ions to small, stationary magnetic island perturbations in a tokamak plasma, as well as the implications for the NTM threshold physics. The effect of finite particle orbit width is substantial. The radial profile of the perturbed ion distribu tion function is shifted radially relative to the magnetic island. This implies that the distribution is no longer flattened across the magnetic island, but instead across a radially shifted drift island. This shift is important for small islands comparable to the ion banana width, in which case a pressure gradient is maintained inside the magnetic island, even if crossfield transport is neglected. The bootstrap current drive for the NTM is then suppressed, with the flows dominated by the electron physics. The response of the electrons to the perturbed elec trostatic potential is such that it tends to heal islands of width w below a critical width w c , thus providing a threshold for NTM growth. We find that, in the absence of other effects, the critical island width scales linearly with ρ θi : w c ∼ 2.76ρ θi .
The new physics of the finite ion orbit width effect is impor tant for a complete theory of the neoclassical tearing mode threshold and, in particular, for designing the NTM control system for ITER. For our theory to fully quantify the NTM theory, we need to address additional physics including the accuracy of the analytic electron response employed here, the finite ion Larmor radius effect, particularly in the vicinity of the island separatrix, as well as the impact of the island rota tion that leads to the ion polarization current. Nevertheless, this work gives a new insight into the physics of small magn etic islands and the NTM threshold.
where, as for the passing particles, α t l is a square matrix and β t l is a vector. The recurrence relations in the trapped region are then: where M (t) l = Q l + P l · α t l+1 . Applying the y boundary conditions at the grid edges allows us to determine α p l=1 and β p l=1 at the deeply passing end (y 1 = 0), and α t l=Ny and β t l=Ny at the deeply trapped end ( y Ny = y max , N y is the number of y grid points). Using the recurrence rela tions (A.4) and (A.6), we can determine all the α l and β l , up to the trapped/passing boundary.
In principle, the solution vectors h p,t l can then be deter mined using equations (A.3) and (A.5), given the solution at the boundary, h l=lc . This can be obtained using the matching conditions [14]: with subscript c corresponding to the grid point at the boundary, and: Here, subscript c + 1 corresponds to the first barely trapped particle grid point ( y c + ∆ y ), c − 1 the first barely passing grid point ( y c − ∆ y ), and so to.
Once h c is determined, then the complete set of solutions in passing and trapped regions can be calculated using the recur rence relations (A.3) and (A.5).

A.1. Momentum conservation and quasi-neutrality
In order to determine the ion response accurately, it is crucial to ensure that the momentum conservation and quasineu trality are imposed selfconsistently. The former is introduced as an additional term in the model collision operator (9): the term in ū (Ḡ i ). This is effectively a weighted parallel velocity moment of the ion distribution function, which adds a degree of nonlinearity. The term in ū i is determined by iterating over the calculation of h i , updating ū i each time until the solution is converged.
Quasineutrality is imposed by determining the perturbed electrostatic potential from ion and electron densities. Given the analytic solution for the electrons [14] summarized in sec tion 3 and the form of the ion distribution function (23), the quasineutrality condition: n i n e implies that: assuming T i = T e . Here, L n0 is the equilibrium density gra dient length scale, δn i is the perturbation in the ion density calculated from h i , and 2π Ω + cos ξ dξ.
The calculation for h i is iterated over until convergence is achieved for Φ, as well as ū i . The converged ion flow is then employed in the electron flow (22). More precisely, the itera tion loop for converging ū i is nested inside the iteration loop for converging Φ. All the steps and flow of the calculations are outlined in the flow chart ( figure A1).