Drift kinetic theory of neoclassical tearing modes in tokamak plasmas: polarisation current and its effect on magnetic island threshold physics

A nonlinear 4-dimensional drift island theory derived in (Imada et al 2019 Nucl. Fusion 59 046016 and references therein) provides qualitative predictions of the plasma response to a stationary neoclassical tearing mode (NTM) magnetic island in a low beta, large aspect ratio tokamak plasma. (Dudkovskaia et al 2021 Plasma Phys. Control. Fusion 63 054001) refines a model for the magnetic drift frequency and exploits the limit of rare collisions, reducing this theory to 3-dimensional and thus providing a more accurate treatment of the trapped-passing boundary layer. The drift island theory is further improved in (Dudkovskaia et al 2023 Nucl. Fusion 63 016020) by introducing plasma shaping and finite beta effects. In the present paper, an improved model is adopted to resolve the drift island separatrix boundary layer, allowing one to investigate the polarisation current contribution that exists around the magnetic island separatrix, including in the presence of the background electric field. In particular, different magnetic topologies from both sides of the separatrix generate a radial discontinuity in the distribution function gradient there, when collisions are neglected. Allowing for collisional dissipation in the leading order distribution function around the separatrix resolves this discontinuity, smoothing the density distribution. The overall effect of the polarisation current on the NTM threshold is then combined from the outer contributions that exist outside the layer, as well as the separatrix layer piece, and self-consistently accounts for the electrostatic potential reconstructed from plasma quasi-neutrality. The corresponding NTM threshold is quantified and compared with previous predictions of (Dudkovskaia et al 2021 Plasma Phys. Control. Fusion 63 054001, Dudkovskaia et al 2023 Nucl. Fusion 63 016020).


Introduction
Neoclassical tearing modes (NTMs) [1][2][3] impose a limit on fusion gain, as well as plasma confinement time in tokamaks.Resulting from filamentation of the current density that flows through the tokamak plasma, they change the equilibrium magnetic topology, leading to a new one, characterised by a chain of magnetic islands.When magnetic islands are large, i.e. much larger than the trapped ion banana orbit width, ρ bi , the plasma thermal pressure is flattened across them.This therefore reduces the total core plasma pressure, resulting in a drop in fusion power [4] in a burning plasma tokamak, like ITER.Furthermore, the pressure profile flattening across the island creates a hole in the bootstrap current density close to the island centre, enhancing the current density filamentation and amplifying the island further.This provides the main driving mechanism for island growth.
The NTM behaviour significantly differs from the above description when magnetic islands are small, i.e. w ∼ ρ bi (a few cm), where w is the magnetic island half-width.Indeed, experimental observations [5][6][7] found that there is some critical (threshold) magnetic island width, 2w c ≈ (2 − 3)ρ bi , below which the pressure gradient is sustained across the island, providing magnetic island 'self-healing'.w c is a key parameter for quantifying the NTM control system (e.g.[8]): the amount of microwave power that must be injected to stabilise the mode, as well as how radially localised it has to be.An accurately quantified w c will help optimise NTM control systems on existing fusion devices, as well as extrapolate these requirements on current drive to design NTM control techniques or avoidance and mitigation strategies for ITER, DEMO and future tokamak power plants.
There are a few approaches aimed to explain the origin of the NTM threshold.Broadly, these can be split between the heat transport threshold model and polarisation-bootstrap threshold model.According to the first one based on [9], the hole in the bootstrap current is the result of the dominant parallel transport along the magnetic field lines, and hence the threshold island width can be estimated by balancing the parallel transport against the perpendicular one, provided the perpendicular gradient length scales are estimated via w.The second model is based on [10] and employs a kinetic theory to calculate the current density perturbation parallel to the magnetic island.In particular, [10] finds w c by balancing the bootstrap current drive and the polarisation current contribution 4 .The polarisation current induced by the magnetic island propagation is found to be stabilising for ω < 0 or ω > ω dia,i (1 + L n /L Ti ) [10], where ω dia,i is the ion diamagnetic frequency and ω is the island propagation frequency with ω < 0 corresponding to the electron diamagnetic direction; L n and L Ti denote the density and ion temperature equilibrium gradient length scales, respectively.The theory of [10] allows for w ∼ ρ ϑi = ε −1/2 ρ bi but requires w ≫ ρ bi , where ρ ϑi is the ion poloidal Larmor radius and ε = r s /R 0 is the tokamak inverse aspect ratio with r s corresponding to the location of the rational surface and R 0 being the tokamak major radius.While [10] quantifies the polarisation contribution that exists outside the magnetic island separatrix, it does not account for a narrow layer in the vicinity of its separatrix [11][12][13][14][15][16][17]  5 .Some previous works (e.g.[13,14]) demonstrated that inclusion of the latter can reverse the effect of the polarisation contribution on the magnetic island stability (which is important for islands rotating with frequencies comparable to the ion/electron diamagnetic frequency) found in [10].In particular, in [13] there is a non-resolved discontinuity in the electron density gradient across the magnetic island separatrix due to different topologies inside and outside the magnetic island which results in a spike in the polarisation current there, large enough to reverse the overall effect of the polarisation current.Inclusion of the cross-field diffusion [17] resolves this discontinuity, smoothing the density distribution across the magnetic island separatrix and decreasing the destabilising contribution to the polarisation current from the separatrix layer.Therefore, whether the polarisation current effect is stabilising or destabilising depends on a subtle balance of these two contributions (one that exists outside the magnetic island and the other one that arises from the thin separatrix layer), and how this balance changes with the island width.
In [18][19][20] we have introduced a drift island formalism to quantify the plasma response to magnetic islands of width close to threshold.According to the drift island formalism, when collisions are rare (with the particle collision frequency being much less than the free streaming along the equilibrium magnetic field lines/characteristic drift frequency), the ions and electrons follow streamlines in phase space.For passing ions and electrons, these streamlines lie on surfaces of the same topology as the flux surfaces of the magnetic island but (1) are shifted radially by an amount that scales with the ion/electron poloidal Larmor radius and, (2) being associated flows perpendicular to the magnetic field (i.e. the polarisation current).Since the polarisation current is not divergence-free, there is also a current that flows parallel to the magnetic field that (1) ensures that the total current is divergence-free and (2) contributes to the magnetic island evolution.This contribution is to be referred to as the 'polarisation' current contribution. 5Note that most of these models consider a simplified (sheared slab) geometry [13,16,17] and/or employ a model electrostatic potential [14,15] of [10].
Table 1.DK (Drift Kinetic)-and RDK-NTM (Reduced Drift Kinetic-NTM) codes and their features.'numerical' denotes the solution of the orbit averaged drift kinetic equation: equation (32) of [19] for DK-NTM, equation (20) of [20] for RDK-NTM (v.1), equation (30) of [21] for RDK-NTM (v.2) and equation (19) of the present paper with included boundary layer effects as discussed below in section 3.1 for RDK-NTM (v.3).'analytic' is based on the analytic results of [10] for the electron response, see [19] for details regarding its numerical implementation in the 4D (DK-NTM) formulation.Spatial variables are pφ, ξ (helical angle) or S; velocity variables are λ (pitch angle), V (particle speed) with σ being the sign of the velocity component parallel to the equilibrium magnetic field.

DK-NTM
RDK-NTM (v. a While both dissipation layers are allowed in the DK-NTM formulation from the mathematical perspective, it is numerically challenging to properly resolve both layers simultaneously when the problem is 4-dimensional. with magnetic drifts, this shift depends on the sign of the parallel velocity component, u, i.e. is in opposite direction for u ≷ 0. We therefore refer to these structures as drift islands.The passing particle distribution function is then found to be a 'drift surface' quantity with its profile being flattened across the drift islands rather than the magnetic island.When w ≫ ρ ϑi , the radial shift of drift islands relative to the magnetic island is negligible, and hence drift island contours almost coincide with the flux contours of the magnetic island, resulting in the density profile flattening across the magnetic island (in accordance with the conventional theory of large magnetic islands, e.g.[2]).At w ≳ ρ ϑi , flattening of the passing ion distribution function across the u ≷ 0 drift islands partially restores the ion density gradient across the magnetic island centre.Note that since ρ ϑe ≪ ρ ϑi ≲ w (ρ ϑe is the electron poloidal Larmor radius), the drift island effect is less significant for electrons, and the electron density profile is still flattened across the magnetic island (in the absence of an electrostatic potential).This difference in the ion and electron density profiles generates a perturbed electrostatic potential to restore the electron density gradient inside the magnetic island, maintaining plasma quasi-neutrality and providing a physics basis for the NTM threshold.Imada et al [18,19] solve a 4D drift kinetic equation for ions (two spatial and two velocity variables) in a low beta, large aspect ratio tokamak plasma, employing the analytic results of [10] for the electron response.The ion and electron responses to the magnetic island are then coupled via the electron-ion collisions and the electrostatic potential required for quasi-neutrality.Reference [20] also exploits the low beta, large aspect ratio tokamak limit, as well as the limit of rare collisions, working in terms of the drift island stream functions and hence solving a 3D drift kinetic equation for ions and electrons (one spatial and two velocity variables) consistent with quasi-neutrality.Note that a transition from the toroidal component of the canonical angular momentum, p φ , employed for the radial variable in [18,19] to the drift island stream function, S, of [20] reduces the problem dimensionality from 4D to 3D for rare collisions.There are certain regions where collisions cannot be treated perturbatively (even in the limit of rare collisions).These are collisional boundary layers around the trapped-passing boundary in velocity space and the drift island separatrix in phase space.The boundary layers are to be introduced as matching conditions, exploiting their thinness (see [10,20] for the trapped-passing boundary layer solution).This 3D formulation with matching boundary layers is numerically beneficial as it allows one to resolve fine structures in boundary layers, not reducing the treatment of the electron response to analytic.The numerical treatment of electrons is numerically more challenging for the 4D model of [18,19] (even when the particle collision frequency is only slightly less than the characteristic drift frequency).To explore these two approaches, two numerical codes have been developed: 4D DK-NTM [18,19] and 3D RDK-NTM (v.1) [20].Their features are compared in table 1.The results of [18,19] and [20] are benchmarked in [20] at plasma collisionalities ν * = 10 −2 and 10 −3 , where the validity regions of the two formulations overlap.
Considering a stationary isolated magnetic island chain in a low beta, large aspect ratio plasma in the absence of the equilibrium radial electric field in the magnetic island rest frame, [18][19][20] focus primarily on bootstrap drive, also capturing polarisation effects due to the perturbed electrostatic potential.The drift island theory is further improved in [21] by incorporating the effects of plasma shaping and finite beta, being referred to as RDK-NTM (v.2) in table 1.In the present paper, we refine the treatment of the drift island separatrix in the 3D formulation of [20,21].Similar to the magnetic island separatrix layer [11][12][13][14][15][16][17], there is a thin boundary layer around the drift island separatrix due to different topologies inside and outside the drift island.When w ≫ ρ ϑi , both magnetic and drift island boundary layers almost coincide.At w ≳ ρ ϑi , the layer around drift islands remains narrow, broadening the layers from both sides of the magnetic island.As was previously noted in [20,21], these separatrix layer effects are expected to be less significant for non-rotating magnetic islands but become crucial when the magnetic island propagation frequency is non-zero (or alternatively in the presence of the equilibrium radial electric field if one works in the magnetic island rest frame).This current version of RDK-NTM is being referred to as RDK-NTM (v.3) in table 1.
The paper is organised as follows.In section 2 we define the magnetic geometry and coordinate system.In section 3 we describe the problem from a mathematical perspective, including the drift kinetic equation for bulk ions and electrons, as well as the corresponding boundary and matching conditions.Its solution then provides the ion and electron distribution functions and the corresponding current density perturbation at each w and ρ ϑi .A series of scans in w at each ρ ϑi then identifies the stationary states of the magnetic island chain, determining the threshold magnetic island width.The corresponding results are presented in section 4 for different tokamak geometries.The impact of the background radial electric field on the polarisation current contribution is also addressed in section 4. The results obtained are summarised and discussed in section 5.

Coordinate system and magnetic topology
We employ the coordinate system of [21] with the equilibrium magnetic field written as where φ is the toroidal angle and ψ is the poloidal flux with the surfaces of constant ψ forming a set of nested toroidal magnetic flux surfaces.I = RB φ with R being the varying tokamak major radius and B φ the toroidal component of the magnetic field.The magnetic field perturbation associated with a single isolated magnetic island chain is introduced as with where b 0 = B 0 /B 0 with B 0 = |B 0 |; ψ = (w ψ 2 /4)(q ′ /q) with w ψ being the island half-width in ψ space related to w in r space via w = w ψ /(RB ϑ ) and B ϑ being the poloidal component of the magnetic field.The helical angle has been defined as ξ = φ − q s ϑ, where q s = m/n is the safety factor, q, evaluated at the rational surface, i.e. the ratio of the poloidal m to the toroidal n mode number; ϑ is the poloidal angle with ∇ψ A prime denotes the derivative with respect to ψ.In [21], it was shown that the second term on the right hand side of equation ( 2) associated with the perpendicular component of the perturbed vector potential, A 1  ⊥ , results in the O δ * ∆ 2 Vt L f 0 corrections in the drift kinetic equation and thus is to be omitted.Here δ * = ρ/L is a small parameter associated with the drift kinetic ordering (ρ is the particle Larmor radius and L characterises the system size) [22], V t is the particle thermal velocity and f 0 is the equilibrium distribution function 6 .Ordering w ∼ ρ bi ≪ r s (ρ bi is the trapped ion banana orbit width and r = r s is the radius of the rational surface associated with the island) ensures the finite orbit width effects of trapped ions are captured, while providing a small expansion parameter, ∆ = w/r s ∼ ρ bi /r s ≪ 1, to simplify the 5D drift kinetic equation (three spatial and two velocity variables).
The spatial coordinates are ψ which labels magnetic flux surfaces, ξ which labels the magnetic field lines on a chosen flux surface and ϑ which measures the distance along a field line.The velocity vector is with The α dependence is removed from the problem due to gyro-averaging in the drift kinetic formalism [22].The velocity coordinates are then represented by the magnetic moment per unit mass, µ, and the total particle energy per unit mass, U: µ = s 2 /2B and U = µB + u 2 /2 + eZΦ /m ≡ V 2 /2 + eZΦ /m.Here eZ and m are the particle charge and mass, V is the particle speed, and Φ is the electrostatic potential.

Plasma response
Reference [21] provides the mathematical details of how we employ the drift kinetic equation to describe the plasma response to the NTM magnetic field perturbation in a realistic tokamak geometry characterised by finite beta values and D-shaped plasmas.A typical drift kinetic problem is 5D, i.e. {ψ , ξ , ϑ, µ, U}.Since we investigate small islands of halfwidth w ∼ ρ bi , the dimension of the problem can be further reduced by employing ∆ as a small parameter to expand the solution (e.g.see [21]).We separate the particle distribution function into the equilibrium Maxwell-Botzmann contribution and a perturbed piece, g = O(∆f 0 ), that describes how plasma responds to the magnetic island: The equilibrium distribution function, f 0 , is assumed to be Maxwellian.The electrostatic potential, Φ, contains the equilibrium piece, Φ eqm = Φ eqm (ψ ), associated with f 0 and perturbation, δΦ = O(∆Φ eqm ), due to the magnetic island: Here ψ = ψ s denotes the location of the rational surface in ψ space, and the subscript 'eqm' indicates the equilibrium quantities.The equilibrium plasma quasi-neutrality is satisfied automatically by the Maxwellian, while δΦ = δΦ (ψ, ξ, ϑ) arises due to the difference in the ion and electron responses to the magnetic island.Noting that the radial variation in equilibrium quantities results in the O δ * ∆ 2 Vt L f 0 corrections in the drift kinetic equation, we evaluate f 0 and the equilibrium ion/electron plasma temperature, T, in equation ( 7) at ψ = ψ s , retaining the radial variation in perturbed quantities only.
Following [18][19][20][21], we work in the island rest frame.Therefore, the Φ ′ eqm (and hence ω E ≡ mΦ ′ eqm /q s ) variation in the island rest frame provides the island propagation frequency, ω, dependence in the reference frame in which the equilibrium electric field is zero far from the magnetic island (see section 2.5 of [23]).This then eliminates the time variation associated with the first order differential operator, i.e. ∂/∂t, from the drift kinetic equation, replacing it with the parametric dependence on ω E .
To leading order in ∆, the toroidal component of the canonical angular momentum, p φ = ψ − ψ s − Iu/ω c with ω c = eZB/m being the cyclotron frequency, is approximately constant on a particle orbit as a consequence of approximate toroidal symmetry.This allows further reduction to 4D, provided p φ is introduced as the radial coordinate instead of the poloidal flux.Specifically, to leading order in ∆, the plasma response is independent of ϑ at fixed p φ .To determine this leading order response, g (0) = g (0) (p φ , ξ, µ, V), we proceed to the following order in ∆, which also contains the higher order plasma response, g (1) , dependent on ϑ7 .To eliminate the g (1) dependent term, we apply the orbit averaging operator, integrating over ϑ at fixed p φ : Here σ ≡ u/|u| = ±1, λ = 2µ/V 2 is the pitch angle and λ c denotes the trapped-passing boundary.λ fin is the upper limit of the λ variation.The bounce points, ϑ = ϑ b , are provided by the condition: u(ϑ b ) = 0.This then provides four orbitaveraged drift kinetic equations for g (0) to be solved in a 4D {p φ , ξ , λ, V} space, i.e. for passing/trapped ions/electrons, respectively, involving (1) different averages employed for passing and trapped particles and (2) different momentumrestoring components in the collision operator for ions and electrons [10].Introducing the dimensionless system of [20,21]: where µ 0 is the magnetic permeability of free space, q ′ s is q ′ evaluated at ψ = ψ s , T e is the equilibrium electron temperature and p eqm is the equilibrium plasma pressure, we write the resulting 4D equation as with the normalised magnetic drift frequency, ωD , being defined as Here β ′ and θ′ have been introduced as β and ϑ differentiated with respect to ψ. Θ denotes the Heaviside step function.
Ĉ on the right hand side of equation ( 10) is the normalised Iu C⟩ pφ ϑ term related to the collision operator, C. Its explicit representation can be found in appendix A or in appendix B of [20] in the low beta limit.Note that in the limit of low beta, large aspect ratio tokamak plasmas, equation (10) with equation ( 11) reduces to equations ( 14), (15)/ equations ( 20), (21) of [20].
Equation (10) is to be solved with the following boundary conditions: (1) at large distances from the magnetic island, the radial derivative of g (0) must match the radial derivative of the equilibrium distribution function, i.e. the island physics, including its impact on the transport processes, is localised around the rational surface; (2) g (0) must be periodic in ξ; (3) to eliminate divergences at the deeply passing and trapped limits in λ, we require a mixed boundary condition there that results from equation (10) (see appendix B).Equation (10) then provides g (0) = g (0) (p φ , ξ , λ, V) = g (0) (ψ , ξ , ϑ, λ, V) for ions and electrons at each ω E .Noting that the electrons have much narrower Larmor orbits (compared to ions) and their motion along magnetic field lines significantly exceeds the guiding centre drifts, the ion and electron responses to the magnetic island perturbation are different.This difference violates plasma quasi-neutrality, requiring an electrostatic potential, δΦ, that must be calculated consistently with the quasi-neutrality condition [10,[18][19][20][21]. 8 This therefore couples electron and ion physics.As discussed in [20], we adopt an iteration process such that at each iteration step in to the electron diamagnetic frequency and x = (ψ − ψs)/w ψ .δn i/e is the perturbed ion/electron density associated with g (0) i/e of equation (7). the perturbed electrostatic potential, δΦ (which is a part of Φ in equation ( 10)) is known.Φ ′ eqm in Φ is an input parameter, and its variation provides the parametric dependence of the solution, g (0) , on the island propagation frequency, ω.
In addition to δΦ and plasma quasi-neutrality, the ion and electron responses are coupled via collisions.Following [10,[18][19][20][21], we employ a collision operator, C, that conserves particles and momentum, equation (62) of [10].This makes equation (10) for ions and electrons integro-differential and couples the electron distribution function to the ion response via the electron-ion collisions.Therefore, we obtain a nonlinear system of 4D integro-differential equations coupled by (1) the momentum conservation term in the electron-ion collision operator and (2) the perturbed electrostatic potential required for quasi-neutrality.
There are certain numerical challenges associated with solving equation (10) in p φ space (being referred to here as the DK-NTM approach) [24], as it relies on the ability of the numerical solver to resolve fine structures in boundary layers.There are two boundary layers: one is in pitch angle space in the vicinity of the boundary between passing and trapped particles, and one that surrounds the drift island separatrix and results from different topologies (and thus averages in helical angle) inside and outside the island.In [20,21] we have introduced a model further reduced in collisions (being referred to here as the RDK-NTM approach): being valid in the limit when the particle collision frequency is much less than the magnetic drift frequency and/or parallel streaming around the island flux contours, it allows one to better resolve boundary layers and thus learn more about the physics related to collisional dissipation, providing an extension to the DK-NTM formulation to low collisionality plasmas (see table 1 for the DK-and RDK-NTM comparison).In this low collision frequency limit, the collision operator can be treated perturbatively, and one then finds that all the spatial variation can be captured in a single stream function, S = S (p φ , ξ, λ, V; σ) so that g (0) = g (0) (p φ , ξ , λ, V) = g (0) (S, λ, V) and the problem reduces to 3D.Introducing collisions perturbatively at the following order and averaging over the drift surfaces, labelled by S, then provides the resulting equation to be solved for g (0)  in S space.Indeed, equation ( 10) is equivalent to where provided the equilibrium quantities are slowly varying functions of ψ that can be evaluated at ψ s and Equation ( 12) can schematically be written as where A denotes the coefficient in front of ∂g (0) /∂ξ S,ϑ,λ,V on the left hand side of equation ( 12) and C represents the right hand side operator of equation ( 12), acting on g (0) .Related to Ĉ, C contains the pitch angle scattering operator, the diffusiontype operator in S space, as well as mixed partial derivative contributions associated with the variable switch from ψ to S [20].In the limit of finite beta, C also includes a contribution that results from the Vlasov operator (second term on the right hand side of equation ( 12)).Following [21], we require where ν * is plasma collisionality and ε is the tokamak inverse aspect ratio evaluated at ψ = ψ s , to allow a perturbative treatment of the right hand side of equation (12).Equation (15) does not allow for any sharp spatial variations in β ′ (consistent with the assumptions made above about the radial variation of equilibrium quantities) and allows for β varying between 10 −2 − 1 for the typical values of ν * in the RDK-NTM formulation, i.e. 10 −4 − 10 −2 .In this limit, to leading order in collisions, confirming that the leading order distribution function, g (0,0) , is independent of ξ at fixed S.
Thus, we learn that the particle distribution function is almost constant along the contours of constant S. Proceeding to next order in collisions, we obtain where g (0,1) is the first order correction in collisions.g (0,1) depends on ξ at fixed S. To eliminate the g (0,1) dependent term, we apply the ξ-averaging operator, ⟨. ..⟩ S ξ , integrating over ξ at fixed S [20], to find i.e. a 3D integro-differential equation that describes how the distribution function g (0,0) varies with S, λ and V at each ω E .Note that equation (18) allows for different collisional models but for the present paper we employ the momentum conserving pitch angle scattering model of [10,[18][19][20][21].While the latter is relatively simple in that it does not include slowing down effects, it is sufficient for the purposes of this paper 9 .Equation ( 18) then produces Here νj is ν ii R 0 /V ti for ions and (νee + ν ei )R 0 /Vte for electrons with V tj being the thermal velocity of species j, and ν jj and ν ei denoting the corresponding collision frequency ( j = i for ions and j = e for electrons).The transport-type coefficients, C SS , C S and C λS (see appendix A for their explicit representation) result from a variable transformation from ψ to S (via pφ to provide additional factorisation, see section D.7 of [23]) and include orbit-averaging at fixed pφ, i.e. they are therefore different for passing and trapped particles in accordance with equation ( 8).The corresponding terms in equation ( 19) describe transport across surfaces of constant S. a(λ, V; σ) and b(λ, V; σ) are the pitch angle scattering weight functions (see appendix A).U ∥j is the momentum-restoring contribution (see appendix A).For ions, U ∥i depends on the perturbed ion distribution function, g (0,0) i , only 10 .For electrons, electron-ion collisions must be retained as well, and hence U ∥e depends on both g (0,0) i and the perturbed electron distribution function, g (0,0) e .Based on equations ( 13) and ( 15), S has different definitions for passing and trapped particles.For example, in the absence of the electrostatic potential, S is proportional to pφ for trapped particles.For passing particles, the surfaces of constant S reproduce the magnetic island flux contours but are radially shifted by an amount ρϑ û + ( Lq/ ŵ)ρ ϑ ωD , proportional to ρ ϑ (see figure 1).Being σ dependent, this shift is in opposite directions for the two streams σ = ±1 (u ≷ 0).Therefore, we refer to these structures defined by the constant S contours as drift islands.Inclusion of Φ calculated from plasma quasi-neutrality is not found to add significant qualitative modifications to the drift island structures [18][19][20][21] but is important quantitatively, and thus its effect is to be retained.Therefore, for trapped and 10 Ion-electron collisions are considered to be small.for passing particles.Here S = Sc corresponds to the drift island separatrix and S min denotes the lower limit of the S variation.At the zeroth iteration in δΦ (i.e.δΦ = 0), Sc and S min can be determined analytically, e.g. they are provided by Sc = ŵ/4 Lq and S min = −ŵ/4 Lq (ω E = 0) and are then to be updated accordingly at the end of each iteration in δΦ, based on the numerical solution for g (0,0) .ξ b1,2 are the maximum and minimum values of ξ on the closed surfaces inside the drift islands (ξ-bounce points) provided by ξ b1,2 = ξ b1,2 S, p φ 0 , λ, V; σ with p φ 0 being a stationary point of S as a function of pφ at each ξ, λ, V and σ. σp φ is then the sign of pφ − p φ 0 .
As was previously shown in [20], the distribution function is flattened across the drift islands (and not the magnetic island).This therefore provides a contribution to the density gradient associated with passing ions at w ∼ ρ ϑi due to the summation over σ in the integral over velocity space.Since ρ ϑe ≪ w, the drift islands for electrons almost coincide with the magnetic island, and the electron density gradient is still flattened across the magnetic island in the absence of a potential.However, this difference in the electron and ion responses drives a perturbed electrostatic potential to maintain plasma quasi-neutrality, restoring the electron density gradient across the magnetic island as well.Note that inclusion of ω E in equation ( 13) impacts the radial shift via Φ ′ eqm ψ =ψ s (ψ − ψ s ) [23], and therefore, depending on its sign, it can further enhance or suppress the drift island effect on the NTM drive.

Boundary layers
The perturbative treatment of collisions becomes invalid in narrow boundary layers: (1) in the region of pitch angle that separates trapped and passing particles (vertical blue area that surrounds λc in figure 2) and (2) in the vicinity of the drift island separatrix (horizontal red area around Sc in figure 2).Due to steep gradients in pitch angle around the trapped-passing boundary, λ = λc, parallel streaming, A ∂/∂ξ | S , there is comparable to collisional dissipation associated with D λ ∂ 2 /∂λ 2 , where D λ is the diffusion-type coefficient in λ space, related to the plasma collision frequency (for definitions, see the following subsection).These steep gradients in λ exist because different orbit-averaging procedures for passing and trapped particles generate a discontinuity at λ = λc when collisions are neglected, which can then be resolved by the collisional layer in pitch angle.The corresponding modifications of the particle distribution function are discussed in [10,20].Another boundary layer arises due to steep gradients in S around the drift island separatrix, S = Sc, where A ∂/∂ξ | S is comparable to the transport term, i.e.D S ∂ 2 /∂S 2 with D S being the diffusion-type coefficient in S space (for definitions, see the following subsection).Indeed, different ξ-averaging procedures at fixed S inside and outside the drift island also generate a discontinuity, and the drift island separatrix boundary layer resolves it.Topologically this is similar to (and replaces) the boundary layers introduced by [9] and [12] for magnetic islands.Indeed, both magnetic and drift island boundary layers almost coincide when magnetic islands are large, i.e. w ≫ ρ ϑi .When w ≳ ρ ϑi , the layer around drift islands remains narrow, therefore broadening the layers from both sides of the magnetic island.Then there is also a special region where both boundary layers overlap, i.e. schematically 2).
3.1.1.λp ⩽ λ ⩽ λc and S t ⩽ S ⩽ Sp region.This is the region where both boundary layers overlap, and the collision operator is dominated by the pitch angle scattering term proportional to S , as well as the diffusion-type term associated with We exploit the narrowness of this region in both S and λ space, provided the plasma has low collisionality, and fix the coefficients in equation ( 22) at λp ≡ λc − ϵ λ [20] and S p/t ≡ Sc ± ϵ S , where λc is the trapped-passing boundary, Sc corresponds to the location of the drift island separatrix, and ϵ λ and ϵ S provide the widths of the trapped-passing boundary layer and the drift island separatrix layer, respectively 11 .Thus, equation (22) reduces to 11 Following [23], and In equations ( 24) and ( 25), the superscript 'in' indicates the phase space region inside the drift island, i.e. evaluated at St, while 'out' corresponds to the region outside the drift island, i.e. evaluated at Sp.Here we have introduced λ = λ − λc and S = S − Sc, which are such that λ = 0 defines the trappedpassing boundary and S = 0 indicates the location of the drift island separatrix.λ < 0 corresponds to the passing region; S ≶ 0 corresponds to the phase space inside/outside the drift island, respectively.Equation ( 23) is to be solved in the region of passing particles for each σ at σp φ = ±1.Equation ( 23) allows for an analytic solution of the Fourier form (to be interpreted similar to equation (97) of [10] where the trapped-passing boundary layer in pitch angle space was investigated): outside and To ensure continuity across S = Sp and S = St, C out 0 and C in 0 are to be provided by drive terms, H + and H − in equation (31) of [20], evaluated at S = S p/t , respectively.The rest of the Fourier coefficients, , are unknown and to be found from matching at S = 0 in inverse k space based on at each σp φ to ensure continuity of g (0) j and its first derivative in S across Sc.Substituting the resulting Fourier coefficients into the above analytic solution, equations ( 26) and ( 27), provides matching across the drift island separatrix, and this is then to be employed to determine the trapped particle solution in the collisional trapped-passing boundary layer, λc ⩽ λ ⩽ λt, as well as to reconstruct the passing solution at λ < λp (see figure 3 for an example layer solution at λ = λp).

λc ⩽ λ ⩽ λ t and S t ⩽ S ≤ Sp region.
Following [10,20], we impose matching conditions provided by equations ( 106)-(108) of [10] at the trapped-passing boundary to ensure continuity of the solution across λ = λc, as well as the same scattering rate of trapped/passing particles into passing/trapped orbits.As explained in [20], S has different definitions from the side of passing (λ < λc) and trapped The equilibrium geometry parameters are chosen as in [20]: inverse aspect ratio at the surface of interest ε = 0.1, plasma elongation κ = 1 and triangularity δ = 0; safety factor at the surface of interest q = 2, normalised safety factor gradient length scale Lq = 1 and the normalised density and ion temperature gradient length scales Ln/ψs = 1 and L Ti /ψs = 1, respectively.Inset: full solution of equation (22).Note that g where neqm is the equilibrium plasma density and V ti is the ion thermal velocity.The red dashed line indicates the location of the drift island separatrix.
According to equations ( 106)-( 108) of [10] and equations ( 26), ( 27) and ( 32), matching at λ = λc provides which then replaces equation (31) of [20] in this region.Equation ( 33) provides a set of three equations for 2N + 1 unknowns, i.e. a t k (N unknowns), b t k (N unknowns) and C t 0 13 .Similar to [20], due to a difference in definitions of ξ and γ t , matching at λ = λc via equation ( 33) cannot be performed in k space.Instead, γ t is a function of ξ, in accordance with equation (31), and hence taking a number of points in ξ space, 3 N, where N must be selected accordingly, i.e.N = 3n, n ∈ N, and treating , providing matching at fixed pφ and ξ.To ensure continuity of the solution in pφ space (i.e. at p φ = p p/t± φ defined in accordance with figure 4), C t 0 is to be provided by the drive term, H t , of equation (31) in [20].Note that H t also ensures that the solution matches the Maxwellian equilibrium far from the magnetic island.Substituting the resulting Fourier coefficients into the above analytic solution in the trapped region, equation (32), provides matching across the trapped-passing boundary, and this is then to be employed as the boundary condition to reconstruct the trapped solution at λ > λt. 13 Similar to drive terms, H ±/t , of [20] and C in/out 0 of section 3.1.1,C t 0 is real.

λp ⩽ λ ⩽ λ t region outside drift island separatrix layer.
In the interval λp ⩽ λ ⩽ λc, S > Sp (outside the drift island) and S min ⩽ S < St (inside the drift island), as well as at the corresponding pφ values from the trapped region, λc < λ ⩽ λt, the collision operator is approximated by a diffusion-type term similar to section 3.1.2and [20]: where D ±/t = (2ν j / V)a λ p/t for passing, σ = ±1, and trapped branches. Here for passing particles outside and inside the Ŝ island.Note that in this region Ŝ = S λ p/t for passing/trapped particles.For trapped particles, γ t is provided by equation (31).
(λ − λc) with the ξ-averaging operator at fixed Ŝ, ⟨. ..⟩Ŝ ξ , being defined as in [20] via an average over ξ for passing particles outside the Ŝ island and trapped particles, and through an integral between ξbounce points, ξ b1,2 , with a weighted sum over σp φ for passing particles inside the Ŝ island.Note that γ ±/t has the same features as the angle variable of [25], providing continuous solutions across ξ = ±ξ b1,2 .The solution of equation ( 34) is of the form: with k , and hence employing again the matching conditions provided by equations ( 106)-(108) of [10], we write which is then equivalent to the system of equation ( 31) of [20].H ±/t is the 0th harmonic, representing the drive terms [20].Due to a difference in the definitions of equations ( 31) and ( 35)-(37), the system of equation ( 39) cannot be solved in k space.Instead, we solve equation ( 39) numerically for a ±/t k , b ±/t k , taking the number of points in ξ space as N ξ = 2N + 1 (N is the number of Fourier harmonics retained in this region) and treating . This provides matching across the trapped-passing boundary at pφ values outside the drift island separatrix (as shown in figure 4).Written in terms of S and evaluated at λ p/t , this solution is then used to recover the passing/trapped particle solution outside the trapped-passing boundary layer (as discussed below in appendix B and in [20]).In the region λp ⩽ λ ⩽ λc, it can be transformed back onto the S grid to provide the C out/in 0 of section 3.1.1 to ensure a continuous solution across S = S p/t .In the region λc < λ ⩽ λt, it provides the C t 0 of section 3.1.2at the values of pφ that correspond to the location of the drift island separatrix layer, pφ = p p/t± φ (see figure 4).

Drift island separatrix layer in the region λ < λp.
In the vicinity of the drift island separatrix outside the collisional trapped-passing boundary layer, St ⩽ S ⩽ Sp and λ < λp, the collision operator is dominated by the diffusion-type term associated with ∂ 2 /∂S 2 λ , and hence which then reduces to Here Din/out S and S are defined in accordance with section 3.1.1.The solution of equation ( 41) is then provided by  42) and (43) to equations ( 26) and ( 27) at λ = λp.Then âin/out k , bin/out k (k > 0) are to be determined from equation (28) at each σp φ to provide continuity of g (0) j and its first derivative in S across Sc.Substituting the resulting coefficients into equations ( 42) and ( 43) provides matching across S = Sc.Equations ( 42) and ( 43) then provide the Dirichlet boundary conditions at S = S p/t to calculate the solution of the drift kinetic equation, equation (19), outside the drift island separatrix boundary layer.

Low collisionality regions
For λ < λp (passing particles, see figure 2) with S > Sp (outside the drift island) and S min ⩽ S < St (inside the drift island), as well as for λ > λt (trapped particles), the collision operator can be introduced perturbatively and the particle distribution function is to be found as the solution of equation (19).Boundary conditions in pitch angle space at λ = 0 and λ = λ fin are mixed and result from equation (19) [23].This ensures that the solution is not divergent at the deeply passing/trapped limit.At the trapped-passing boundary, λ = λc, we match the passing and trapped branches of the distribution function in accordance with section 3.1.At large distances from the magnetic island, the radial derivative of the solution must match the radial derivative of the equilibrium distribution function.This is provided by the Neumann boundary condition in S space [23] at S → +∞ (or pφ → ±∞, see figures 3 and 4), which is typically a few w away from the magnetic island.To ensure that the solution is continuous in pφ space across the drift island centre (i.e. the drift island O-point, pφ = p φ 0 at each λ, V and σ, see figure 4 for example), we match the σp φ = ±1 branches of the passing particle distribution function (see figure 3) in S space at S = S min .At the drift island separatrix, S = Sc, we match the branches of the passing particle distribution function that exist inside and outside the drift island at each σ in accordance with section 3.1.
We impose the Dirichlet boundary condition at S = Sp and S = St to ensure the passing particle solution is continuous across the drift island separatrix layer.In pitch angle space, similar to [20], we impose the Dirichlet boundary condition at λ = λp and λ = λt, provided by the boundary layer matching solution of section 3.1.Discretising g (0,0) of equation ( 19 (44) at each λ grid point, j, and calculating α σ,p/|σ|,t j and β σ,p/|σ|,t j at j = 0 (passing particles, λ = 0) and j = N p2 (trapped particles, λ = λ fin , see figure 2) from the boundary conditions at the deeply passing/trapped limit 14 , we find all α σ,p/|σ|,t j s and β σ,p/|σ|,t j s at each j up to λ = λp from the side of passing particles and back to λ = λt for trapped particles.Substituting the boundary layer solution obtained in section 3.1 and evaluating it at λ = λ p/t in equation ( 44) then provide g σ,p/|σ|,t j at each λ grid point between 0 ⩽ λ < λp and λt < λ ⩽ λ fin .The details of the numerical solution are provided in appendix B. Note that first we iterate over U ∥j of equation ( 19) to converge the momentum-restoring contribution in the collision operator.This is inside the second iteration loop to converge the perturbed electrostatic potential.
An example solution is shown in figure 3 as a function of S and in figure 4 as a function of pφ.For comparative purposes, in appendix C we benchmark the solution of equation ( 19) with boundary layer effects included against (1) the DK-NTM solution of [20] and (2) the RDK-NTM solution of [20] which is obtained in the absence of the drift island separatrix layer physics but captures the trapped-passing boundary layer effects.In particular, in appendix C we consider the ion 'density moments' defined via and the 'flow moments' via Note that the equilibrium plasma parameters in appendix C are chosen as in [20], i.e. assuming a low beta, large aspect ratio tokamak plasma.

Parallel flows in the presence of the magnetic island
While inclusion of the drift island separatrix layer physics is not found to significantly modify the density moments (see appendix C), it qualitatively changes the shape of the flow moments (and hence the flows parallel to the magnetic field) in the (ψ, ξ) plane (e.g.see figures 5 and 6), compared to the results of [20] even in the absence of the equilibrium radial electric field (ω E = 0).Indeed, as the island width decreases towards ρ ϑi , we previously observed some negative (stabilising) flows that spread into the magnetic island region near the X-points (ξ = ±π in figure 6(a)) [20].These flows are then found to compete with growing positive (destabilising) flows outside the magnetic island separatrix, contributing to the NTM threshold.From figure 6(a), it is apparent that these flows are due to part of the parallel flow that varies on the flux surface and flux surface averages to zero which is typically associated with the ion 'polarisation' flow.Since ω E = 0 in this case, the only polarisation effect that is present here is due to the perturbed electrostatic potential.The latter is expected to be sensitive to how the solution is treated in the   vicinity of the island separatrix.In particular, in [20] we considered the diffusion-type terms in S space as contributions ∼ νj times smaller than free streaming (i.e.∆ ≲ νj ≪ 1) in the drift island separatrix layer, which is refined in the present paper in accordance with section 3.1, updating the result for the ion parallel flow (figure 6(b)).
In figure 6(b), we see that the shape of the ion parallel flow contours is modified by the more accurate treatment of the drift island boundary layers: (1) negative (stabilising) flows surrounding X-points and spreading into the magnetic island in figure 6(a) are replaced with negative (stabilising) flows inside the magnetic island separatrix in figure 6(b), in addition to (2) the positive (destabilising) flows outside the separatrix being further amplified compared to figure 6(a).Therefore, in figure 6(b) we observe a competition between the outer positive flows and inner negative flows that determines the NTM threshold.This is then similar qualitatively to the behaviour observed in [20] (figure C8(b)) for the DK-NTM ion flow contour at the plasma collisionality, ν i = 10 −3 .
While inclusion of the drift island separatrix layer physics changes the shape of the flow moments and flows parallel to the magnetic field line, in the following section we demonstrate that it has a relatively small impact on the NTM threshold predictions made in [20] based on the RDK-NTM approach for a low beta, large aspect ratio tokamak and in [21] with plasma shaping and finite beta effects being incorporated.

Neoclassical drive. Polarisation contribution
In the previous section we described the procedure to calculate the leading order ion and electron responses to a single isolated chain of magnetic islands, as well as the corresponding perturbed electrostatic potential required for quasineutrality.Velocity space integration with appropriate weight functions then provides the current density perturbation parallel to the magnetic field, j ∥ , based on these ion and electron responses.The magnetic island growth/ decay rate is then reconstructed from Ampère's law [10].Schematically, this reads (2τ R /r 2 s )dw/dt = ∆ ′ (w) + ∆neo (w) [20], where τ R is the resistive diffusion time with ∆ ′ being the classical tearing mode stability parameter determined by the global equilibrium properties [26] and ∆neo incorporating all the local neoclassical tokamak physics, which primarily includes bootstrap, polarisation and curvature effects [27,28]: where j∥ is the ϑ average of the current density perturbation, j ∥ .Note that generally there might be additional contributions to the magnetic island evolution equation, such as those associated with the NTM control systems [8], coupling to other modes [29][30][31] or additional nonlinear corrections to ∆ ′ associated with the equilibrium current density profile [32,33] when magnetic islands are much larger than the resistive layer width (not relevant to the situation considered in the present paper).The NTM threshold island half-width is then obtained as a stationary point of dw/dt as a function of w, i.e. when ∆ ′ is balanced against ∆neo: ∆ ′ + ∆neo( j ∥ ) = 0. Defining the 'polarisation' current density as part of j ∥ that varies on the magnetic island flux surface (and averages to zero over a flux surface) and 'bootstrap' current density as part of j ∥ that is constant on a flux surface, we split ∆neo = ∆ pol + ∆ * bs with respectively.Here ⟨. ..⟩Ω ξ represents the magnetic island flux surface average (e.g.see equation (35) of [20]): with Ω = 2(ψ − ψ s ) 2 /w 2 ψ − cos(nξ) denoting the magnetic island flux contours, in agreement with equations ( 2) and ( 3).Note that ∆ * bs generally includes additional effects, such as effects of plasma shaping, and reduces to the conventional bootstrap drive term proportional to 1/w in the limit of large w (e.g.[34]).With the assumptions made in [20,21], ∆ * bs reduces to a combination of bootstrap and curvature contributions, ∆ bs + ∆cur, of [20,21].
Note that the drift kinetic theory of [18][19][20][21], including the present paper, quantifies ∆neo only.∆ ′ being governed by the global equilibrium properties is typically derived from ideal MHD and hence can be considered as an input parameter in our formalism 15 .Therefore, instead of evaluating the experimental marginal island width, we introduce the threshold as a critical island half-width, wc, below which the localised NTM drives provided by ∆neo are stabilising.We note though that ∆neo is typically much larger than ∆ ′ at island widths that differ only slightly from 2wc, which then almost removes the distinction between 2wc calculated here and the actual experimental threshold (e.g.see [7]).
In figure 8 we plot ∆ pol as a function of ω E , which characterises the equilibrium radial electric field.Note that this ω E dependence also describes how ∆ pol depends on the island propagation frequency, ω, in the reference frame where the equilibrium radial electric field is zero at large distances from the magnetic island.Figure 8 is qualitatively close to figure 2 of [14] that shows the actual (perpendicular) polarisation current around the magnetic island as a function of ω (see footnote 4).First, ∆ pol = ∆ pol (ω) is approximated by a parabolic function at large ω, which was previously reported, e.g. in [10] (with ∆ pol being stabilising in the absence of the magnetic island separatrix layer physics) and [14] (with ∆ pol being destabilising with the separatrix layer incorporated).Second, ∆ pol is a linear function of ω around ω = 0.In [14,15] the latter, as well as the consequent sign reversal of the polarisation current, is explained by the competition between the toroidal drift of trapped ions and electric drift.However, the sign reversal at ω ≈ ω dia,e is more interesting (ω E = −0.89ωdia,e in figure 8) as it is in a frequency range characteristic of island rotation (e.g.[10,36]).The value of ω E , at which this sign reversal occurs, is not found to change significantly with w/ρ ϑi .Therefore, whether the polarisation contribution is stabilising or destabilising is sensitive to the precise value of the island propagation frequency.
In figure 9 we plot ∆neo, ∆ * bs and ∆ pol as functions of w for different values of the ion poloidal Larmor radius.Figure 9(a) is based on the RDK-NTM model of [20], while figure 9(b) captures the drift island separatrix layer physics in accordance with section 3.1.Strictly speaking, [20] considered the diffusion-type terms associated with ∂ 2 /∂S 2 λ as corrections of order νj ≪ 1 for the entire range of S variation, including the drift island separatrix layer, therefore capturing the destabilising effects of the polarisation currents, but only partially.As one would expect from previous theory models (e.g.[14]) and as can be seen from figure 9, allowing for the S diffusiontype terms being comparable to free streaming enhances the destabilising influence of ∆ pol near the threshold.The latter changes the balance between ∆ pol and ∆ * bs , slightly decreasing the NTM threshold in figure 9: from 0.78ρ ϑi (figure 9  Ratio of the critical island width, 2wc, to ρ bi plotted as a function of plasma triangularity, δ, for the equilibrium parameters employed in [21] (figure 5(d)) with ε ≈ 0.1 and plasma elongation κ ≈ 1 at the q = 2 rational surface.'RDK-NTM' corresponds to the result of figure 5 of [21] obtained in the absence of the drift island separatrix layer effects.'RDK-NTM (layer)' indicates the RDK-NTM result obtained in the present paper with included boundary layer effects.Adapted from [21].© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.
of the separatrix layer physics changes the balance between stabilising and destabilising flows (compare figures 6(a) and (b)), and hence the balance between ∆ * bs and ∆ pol (see figure 9), slightly decreasing the NTM threshold in figure 9. Positive triangularity amplifies further this destabilising effect of the separatrix layer close to threshold (e.g.see figure 10 where we show an artificial δ parameter scan; note that here δ is the triangularity of the flux surface where the NTM is located, and the value of δ consistent with the equilibrium parameters at the surface of interest is δ ≈ 0.1 in figure 10).
In figure 11 we show the ratio of the critical magnetic island width to ε 1/2 ρ ϑi at certain values of ε 1/2 .Experimental 16 Reproduced from LA HAYE R.J., BUTTERY R.J., GERHARDT S.P., SABBAGH S.A. AND BRENNAN D.P. 2012 ASPECT RATIO EFFECTS ON NEOCLASSICAL TEARING MODES FROM COMPARISON BETWEEN DIIII-D AND NATIONAL SPHERICAL TORUS EXPERIMENT PHYS.PLASMAS 19 062506, https://aip.scitation.org/doi/10.1063/1.4729658;FIGURE 6, with the permission of AIP Publishing.The figure has been adapted from the original.Table 2. Some equilibrium quantities used in figures 11 and 12.Note that the elongation shear, sκ, and triangularity shear, s δ , are defined as in [37].All equilibrium quantities are evaluated at the flux surface of interest, q = 2, in accordance with figure 12  data reproduces figure 6 of [7], where effects of finite aspect ratio were studied at DIII-D and NSTX.The rest is the theory prediction based on the RDK-NTM formalism (1) without plasma shaping and separatrix layer effects based on [20], (2) extended to a finite beta, finite aspect ratio tokamak limit of [21] and (3) described in the present paper with plasma shaping included, as well as the boundary layer physics.The ⋆ marker denotes the RDK-NTM scaling of [20] obtained in the limit of small inverse aspect ratio.The −⋆ and ⋆− markers indicate the RDK-NTM prediction for the geometry parameters characteristic of [7]: elongation κ ≈ 2, inverse aspect ratio ε ≈ 0.3 without (δ = 0) and with (δ ≈ 0.4) triangularity effects, in close agreement with experiment.In figure 11 we also show the RDK-NTM predictions for three different aspect ratios to investigate the impact of plasma shaping, as well as the separatrix layer effects, in a wider range of ε variation.The latter is supported by the equilibrium reconstructions of figure 12 and table 2. We note that the corresponding KSTAR and MAST-U experiments were not designed to measure the threshold island width, and hence figure 11 shows only the theory prediction for these equilibria.According to figure 11, the difference between (1), ( 2) and ( 3) is minimal when the inverse aspect ratio (and also triangularity) is small and grows with increasing ε and δ.Therefore, we have not found any significant quantitative difference in the wc variation with poloidal beta obtained in [21] for the EAST plasma with ε = 0.13 at the q = 2 rational surface (e.g.see figure 13).
Inclusion of finite aspect ratio effects visibly decreases the NTM threshold prediction at ε > 0.3, in agreement with [21].Interestingly, a resolved separatrix layer decreases the threshold even further at larger ε (e.g.see ε 1/2 ≈ 0.7 in figure 11). 17Experimental magnetic island width threshold is reconstructed as αψ schematically.Here α is a constant value provided by ECE (electron cyclotron emission) [41,42], and hence it accounts for the distance between the two ECE channels, which is less than 2cm in EAST experiments, resulting in a constant error value of less than 40%.Being constant, α does not impact the qualitative behaviour of 2wc as a function of β ϑ observed in experiment but might quantitatively influence the entire functional dependence.ψr, in contrast, results from the noise level of diagnostics (saddle coils), and hence is responsible for error bars (random error) of a few % in the experimental data of figure 13.

Discussion and future work
The theory and simulations presented here have primarily been motivated by the importance of the separatrix layer physics [11][12][13][14][15][16][17] in quantifying effects associated with polarisation current and hence determining the threshold island width.We have therefore extended the RDK-NTM approach of [20,21] to incorporate the physics associated with a narrow collisional boundary layer that exists around the drift island separatrix and resolves a discontinuity in the drift island stream function S derivatives that arises due to different topologies (and hence ξ-averaging) inside and outside the drift island.
As explained earlier, the DK-NTM approach of [18,19] captures both the trapped-passing boundary layer, as well as the separatrix boundary layer, from a mathematical perspective.However, being a 4D problem, it depends on the ability of the DK-NTM code [20] to resolve these boundary layers, which is challenging numerically and therefore also requires the analytic treatment of the electron plasma component.Introducing the drift island formalism [20] and working in terms of the drift island stream functions, S, when collisions in a plasma are rare (plasma collisionality ν * ≲ 10 −2 ), reduces the problem dimensionality to 3D, which then allows a more accurate resolution of the boundary layers, as well as the numerical drift kinetic treatment of the electron distribution function.
References [18][19][20] quantify the bootstrap drive, also (partially) capturing polarisation effects that arise due to the perturbed electrostatic potential.Dudkovskaia et al [21] extends [18][19][20], introducing plasma shaping and finite beta effects.In addition to already quantified bootstrap and curvature contributions to island evolution, the present paper captures polarisation effects: those that exist outside the magnetic island separatrix, as well as the polarisation contribution that results from the separatrix boundary layer.This therefore provides the final ingredient in quantifying ∆neo = ∆ * bs + ∆ pol , i.e. a contribution to the magnetic island evolution that is associated with the local resistive layer currents 18 .
Based on the results presented, we summarise • Inclusion of the drift island separatrix layer physics qualitatively changes the shape of the passing ion flow moment contours and hence contours of the ion parallel flow, compared to the RDK-NTM results of [20].In particular, it changes part of the parallel flow that varies on the magnetic island flux surface, i.e. the 'polarisation' flow.The latter was anticipated (see section 7 of [20]) since the polarisation flow is known to be sensitive to the separatrix layer physics and the perturbed electrostatic potential there (e.g.compare [13] with a non-resolved discontinuity in the electron density gradient and [17] where the discontinuity is resolved by the cross-field diffusion).
• This is then found to change the balance between the part of the current density perturbation parallel to the magnetic field that varies on the magnetic island flux surface (and hence ∆ pol ) and the part that is constant on a flux surface (and hence ∆ * bs ), modifying the total contribution associated with the resistive (neoclassical) layer currents, ∆neo, compared to the RDK-NTM results of [20,21].
• When effects of the magnetic island separatrix are included, the polarisation contribution, ∆ pol , is found to have a relatively complex dependence on the background radial electric field (or alternatively on the island propagation frequency in the reference frame where the equilibrium electric field is zero), e.g. in agreement with [14].It is important to note that ∆ pol reverses its sign at frequencies comparable to the electron diamagnetic frequency, i.e. the impact of ∆ pol on magnetic island stability depends on the exact value of the island propagation frequency.• At low inverse aspect ratio (ε < 0.3) (characteristic of conventional tokamaks) and low triangularity, the impact of the drift island separatrix layer physics on ∆neo and thus on the threshold magnetic island width is found to be minimal.Therefore, the RDK-NTM prediction of [20], i.e. 2wc = 2.85ρ bi , is reliable in this region (e.g.see figure 11 for the EAST equilibrium reconstruction shown in figure 12).• As we approach ε ≈ 0.3 (figure 11) and/or increase δ (figure 10), effects of plasma shaping start playing a role, decreasing the radial shift of the drift islands [21] and reducing the threshold.Inclusion of the separatrix layer effects broadens and amplifies the destabilising ion flows outside the magnetic island separatrix (figure 6), reducing the threshold further (e.g.see figure 11 for the MAST-U equilibrium reconstruction shown in figure 12).
The (R)DK-NTM theory provides ∆neo = ∆neo(w), i.e. the contribution to magnetic island growth/ decay due to neoclassical currents localised to the vicinity of the rational surface.This can then be used directly as an input to quantify the NTM control system, in particular the NTM threshold island width in terms of the electron cyclotron deposition width (for the electron cyclotron NTM stabilisation) to calculate the NTM stabilisation criteria.An alternative approach is to include the current drive source term in the electron drift kinetic equation in addition to the momentum conserving collision operator of [10] to enable plasma-wave interactions [43].
The theory presented here is subject to certain limitations.First, it describes bulk, thermal ions and electrons and how their distribution functions are influenced by the presence of the magnetic island.To make the theory representative of a burning tokamak plasma (e.g. as in ITER or STEP), one needs to introduce a population of fusion produced alpha particles.Second, the drift kinetic approach generally limits the applicability of the (R)DK-NTM formalism, imposing a limit on B ϑ /B 0 (as discussed in [21]).The latter makes any drift kinetic theory potentially questionable for a spherical tokamak plasma characterised by B ϑ ≲ B 0 , including the 'MAST-U equilibrium' case of figure 11 (see table 2 for B ϑ /B 0 ).To allow for B ϑ ∼ B 0 , one needs to retain the finite ion Larmor orbit effects 19 , which therefore requires the theory extension to gyrokinetics.Extension to gyrokinetics is also desirable as finite Larmor radius effects might potentially be important in the separatrix layer, and it will allow one to investigate the impact of plasma turbulence on the NTM physics close to threshold when the time dependent small scale fluctuations are included.Third, equilibrium quantities have been considered slowly varying functions of ψ, which might be violated in certain plasma regions, particularly in spherical tokamaks.To allow for sharper equilibrium radial profile variations and associated equilibrium neoclassical currents that exist in the absence of the NTM magnetic island, one would need to retain departures from a Maxwellian in the equilibrium distribution function up to O(δ 2 * f 0 ) (see section 4 of [44]), as well as fully capture contributions due to ∇ ψ × b 0 in the NTM magnetic field perturbation (e.g.see [45] for the radial variation of the perturbed flux function) together with resolving B ϑ ∼ B 0 for consistent ordering.The ∇ ψ × b 0 contribution in the later case will also likely require a more realistic (radially asymmetric [46,47]) magnetic island geometry.
European Commission.Neither the European Union nor the European Commission can be held responsible for them.

Appendix A. Right hand side operator in S space
Following [20] and employing the collision operator of [10], we find where with R = R/R 0 , where R 0 is defined as in [37,48].Note that λB 0 is dimensionless.The momentum-restoring contribution, U ∥j , is defined as for ions and for electrons.Here Vj = V j /V tj is the normalised ion (j = i)/ electron (j = e) velocity, and νee = νeeR 0 /Vte and νei = ν ei R 0 /Vte with νe = νee + νei .In the low beta limit, equation (A.1) provides the results of appendix B of [20].

Appendix B. Numerical algorithm in low collisionality regions
Employing the finite difference scheme for the derivatives with respect to λ [23], we rewrite equation ( 19) in a matrix form: Note that y = (S − S min ) 1/2 for passing particles and y = S for trapped particles in accordance with equation (13).The matrix elements of P σ,p/|σ|,t j , Q σ,p/|σ|,t j and R σ,p/|σ|,t j contain the appropriate coefficients of equation (19) and boundary conditions in S (or y) space (see [23]), as well as spacing in λ and y direction.

B.2. Mixed boundary condition at λ = λ fin (deeply trapped limit)
Similar to λ = 0, for deeply trapped particles at λ = λ fin ( j = Np 2 , see figure 2 Np2 ) arises due to the momentum conservation term in the collision operator.In contrast to passing particles, the contribution of trapped particles to the momentum conservation integral is small primarily because of the sum over σ in the orbit averaging operator at fixed pφ for trapped particles, i.e. their distribution function g (0) j (pφ , ξ , λ, V) = g (0) j (ψ , ξ , ϑ, λ, V; σ) with pφ = pφ (ψ , ϑ, λ, V; σ).Then the layer solution in the vicinity of the trappedpassing boundary, equation (38) with the Fourier coefficients provided by equation ( 39), evaluated at λ = λp (for S min ⩽ S ⩽ St and S > Sp) together with α σ,p j and β σ,p j at each j up to Np 1 (see figure 2) provide g σ,p j at each j back to j = 0 via equation (B.4).For trapped particles, equation (38) evaluated at λ = λt together with α |σ|,t In figure C1 we plot the DK-and RDK-NTM ion density moments as functions of pφ for two different ξ values in the passing region of λ variation for the equilibrium parameters of [20].The DK-NTM results are those presented in [20] for a low beta, large aspect ratio tokamak plasma.In figure C1 we also show two RDK-NTM solutions: obtained in [20] with non-resolved separatrix layer and obtained in the present paper that accounts for the boundary layer physics.As can be seen from figure C1, the impact the separatrix layer has on the ion density moments is not particularly large but it becomes crucial for flow moments, especially near the magnetic island separatrix well in the passing region (figure C2) and at w/ρ ϑi ≳ 1 (figure C3).

Similarly, employing
At w/ρ ϑi ≫ 1 (figure C2) the ion flow moment is expected to be zero across the magnetic island, in accordance with the conventional theory of large magnetic islands with no sinks or sources within the island.The DK-NTM result of [20] used here for benchmarking slightly differs from this conventional prediction at ρ ϑi ≪ w in the deeply passing limit (figures C2(a) and (d)), which can be anticipated due to larger numerical errors in the 4D DK-NTM model at this length scale (as discussed in appendix C of [20]).Apart from the deeply passing limit (λ = 0), the flow moments presented in figure C2 are in close agreement.Note that inclusion of the separatrix layer effects significantly reduces the jump in the RDK-NTM ion flow moment across the separatrix (figures C2(a)-(c)) compared to the case of [20] (figures C2(d)-(f )).
In figure C3 we plot the ion flow moments as functions of ψ at two different ξ values, while varying λ in the passing region.The RDK-NTM result here includes the resolved separatrix boundary layer.The flow moments in figure C3 agree well (figures C3(a), (b), (d) and (e)) up to the trapped-passing boundary layer (figures C3(c) and (f )).While flow moments in figures C3(c) and (f ) still agree reasonably well, there are certain small discrepancies in the DK-and RDK-NTM results, in particular in the vicinity of the magnetic island separatrix (ψ ≈ ±0.02 in figure C3).As discussed in [20], it is computationally challenging to resolve both layers (one that surrounds the trapped-passing boundary and one in the vicinity of the island separatrix) simultaneously in DK-NTM with the 4-dimentional solution.Note that ξ = 0/ ξ = ±π corresponds to the magnetic island O-point / X-point.'RDK-NTM' denotes the RDK-NTM solution of [20] in the absence of the drift island separatrix layer effects.[20] in the absence of the drift island separatrix layer effects.'RDK-NTM (layer)' captures the drift island separatrix layer via matching described in section 3.1.

Figure 1 .
Figure 1.Magnetic island flux contours (in the middle) and contours of constant S (i.e.drift island flux contours) for σ = +1 (right) and σ = −1 (left).The green contour shows the separatrix of the magnetic island.The red/blue dashed contour shows the σ = ±1 drift island separatrix, shifted by a value proportional to ρ ϑ relative to the magnetic island.

Figure 2 .
Figure 2. Schematic representation of the solution technique (see appendix B).

Figure 4 .
Figure 4. Leading order ion distribution function, g (0,0) i , of figure 3 plotted as a function of pφ.Inset: zoom in the vicinity of the drift island separatrix based on the full solution of equation (22).The green dashed line indicates the equilibrium distribution function in the absence of the magnetic island perturbation.Note that for passing particles (σ = ±1), p p/t+ φ corresponds to S p/t at σp φ > 0, while p p/t− φ corresponds to S p/t at σp φ < 0.
0).To ensure continuity of the solution across λ = λp, Ĉout 0 and Ĉin 0 are to be provided by C out 0 and C in 0 of section 3.1.1,respectively, via matching equations (

Figure 6 .
Figure 6.Colour contours of the ion parallel flow, u ∥i (10 −2 ), in the (ψ, ξ) plane at ϑ = −π; w/ρ ϑi = 2.86, ν * i = 10 −3 (plasma collisionality).u ∥i is normalised to the ion thermal velocity.The equilibrium geometry parameters are chosen as in [20]: inverse aspect ratio at the surface of interest ε = 0.1, plasma elongation κ = 1 and triangularity δ = 0; safety factor at the surface of interest q = 2, normalised safety factor gradient length scale Lq = 1 and the normalised density and ion temperature gradient length scales Ln/ψs = 1 and L Ti /ψs = 1, respectively.The green curve indicates the position of the magnetic island separatrix.The corresponding perturbed radial electric field is shown in figure 7. (a) Adapted from [20].© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.

Figure 7 .
Figure 7. Colour contours of the normalised perturbed electrostatic potential, δΦ (10 −3 ), differentiated with respect to ψ and plotted in the (ψ, ξ) plane for the parameters of figure 6. δΦ is normalised to Te/e.Note that here ω E = 0.The green curve indicates the position of the magnetic island separatrix.

Figure 9 .
Figure 9. Neoclassical contributions, ∆neo = ∆ * bs + ∆ pol , plotted as a function of w at different ρ ϑi ; ω E = 0, ν i = 10 −3 (plasma collisionality), ε ≈ 0.1 based on (a) the RDK-NTM model of [20] in the absence of the drift island separatrix layer effects and (b) the RDK-NTM result obtained in the present paper with included boundary layer effects.Red, green and purple curves correspond to ∆neo, ∆ * bs and ∆ pol , respectively, with circle, square and diamond markers corresponding to ρ ϑi = 10 −3 , 2 × 10 −3 and 3 × 10 −3 .The cyan '+' curve in (b) indicates the prediction of (a) for ∆neo.The blue dashed/solid curve in (a) indicates the analytic theory prediction of large islands for the bootstrap/polarisation contribution.Note that ∆neo < 0 is stabilising; w, ρ ϑi are normalised to rs.(a) Adapted from [20].© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.

Figure 10 .
Figure 10.Ratio of the critical island width, 2wc, to ρ bi plotted as a function of plasma triangularity, δ, for the equilibrium parameters employed in [21] (figure 5(d)) with ε ≈ 0.1 and plasma elongation κ ≈ 1 at the q = 2 rational surface.'RDK-NTM' corresponds to the result of figure 5 of [21] obtained in the absence of the drift island separatrix layer effects.'RDK-NTM (layer)' indicates the RDK-NTM result obtained in the present paper with included boundary layer effects.Adapted from [21].© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.

Figure 12 .
Figure 12.Magnetic flux surface contours for the equilibrium used in figure11.The experimental equilibrium is fitted to the Miller parametrisation[37].For a detailed description of the Miller parametrisation process, see appendix C of[38].The EAST experimental equilibrium reconstruction is based on the EAST discharge number 91972 (5300 ms)[21,39].The KSTAR experimental equilibrium reconstruction is based on the KSTAR discharge number 32129 (2300 ms).The MAST-U experimental equilibrium reconstruction is based on the MAST-U discharge number 45272 (475 ms)[40].
to the matrices in equation (B.1), containing the appropriate coefficients of the drift kinetic equation and grid spacing in λ and y. the momentum conservation term in the collision operator.Note that we use second order central differencing formulas for the equation and second order one-sided differencing formulas for the boundary conditions (backward for passing branch and forward for trapped branch).Assuming a linear relation between g σ,p j at jth and ( j + 1)th grid points, we write g of passing particles, λ ⩽ λp.Here α σ,p j is a square matrix of size Ny × Ny and β σ,p j is a vector of length Ny.Combining equations (B.1) and (B.4), we obtain the following recurrence relation: α λt ⩽ λ ⩽ λ fin provide g |σ|,t j up to j = Np 2 via equation (B.10).Appendix C. Benchmarking of the (4D) DK-and (3D) RDK-NTM (with/without drift island separatrix layer matching) approaches

Figure C2 . 2 (
Figure C2.The ion flow moments, comparing the results of (R)DK-NTM models, plotted as a function of ψ for different λ (similar to figure C1) across the magnetic island O-point, ξ = 0, at
Ny in the passing/trapped region at each λ grid point, Ny associated with the momentum conservation term in the collision operator.Ny is the number of points in the y direction: Ny = N in y at S min ⩽ S < St (λ < λp), 1)for passing particles (σ = ±1) inside (S min ⩽ S < St) and outside (S > Sp) the drift island and y at S > Sp (λ < λp) and Ny = N t y at λ > λt.
), the boundary condition is given by Np2 are defined similarly to the matrices in equation (B.2), containing the appropriate coefficients of the drift kinetic equation, boundary conditions in S (or y) space and grid spacing in λ and y.