a repository of Nonlinear Kinetic Ion Response to Small Scale Magnetic Islands in Tokamak Plasmas : Neoclassical Tearing Mode Threshold Physics

A new drift-kinetic theory of the ion response to magnetic islands in tokamak plasmas is presented. Small islands are considered, with widths w much smaller than the plasma radius r , but comparable to the trapped ion orbit width ρ bi . An expansion in w=r reduces the system dimensions from five down to four. In the absence of an electrostatic potential, the ions follow stream lines that map out a drift-island structure that is identical to the magnetic island, but shifted by an amount ∼ few ρ bi . The ion distribution function is flattened across these drift islands, not the magnetic island. For small islands, w ∼ ρ bi , the shifted drift islands result in a pressure gradient being maintained across the magnetic island, explaining previous simulation results [E. Poli et al. , Phys. Rev. Lett. 88 , 075001 (2002)]. To maintain quasineutrality an electrostatic potential forms, which then supports a pressure gradient in the electrons also. This influence on the electron physics is shown to stabilize small magnetic islands of width a few ion banana widths, providing a new threshold mechanism for neoclassical tearing modes — a key result for the performance of future tokamaks, including ITER.

Magnetized plasmas are susceptible to tearing mode instabilities. These are characterized by the evolution of magnetic islands, which arise from a filamentation of the component of current density along magnetic field lines. The change in magnetic topology associated with these islands has an impact on the confinement of the plasma by the magnetic field. It is therefore important to determine the conditions under which they grow to large amplitude. To address this, it is necessary to understand how ions and electrons respond to magnetic islands, what currents that response creates, and whether those currents act to amplify or heal the island.
In the simplest picture, particles free-stream along magnetic field lines. As a result, their distribution functions are constant on the perturbed magnetic flux surfaces of the island. In the absence of heat and/or particle sources, this results in a flattening of the distribution function across the island O point. In this Letter, we show that the particle drifts have a significant impact on this picture, especially when the width of particle orbits associated with those drifts are comparable to the island width. As a particular example we focus on the tokamak, which provides a good illustration of the effect because (a) it has E × B, grad-B, and curvature drifts, and (b) the results have consequences for an important tokamak instability called the neoclassical tearing mode (NTM).
In tokamak plasmas, the current density filamentation that drives NTMs is typically dominated by the bootstrap current. In toroidal geometry, a fraction of particles are trapped in the region of low magnetic field, executing closed, banana-shaped orbits. With the presence of a pressure gradient, the finite banana width of those trapped particle orbits drives opposing flows in the ions and electrons along magnetic field lines, in a similar mechanism to that responsible for diamagnetic flows. This seeds the aforementioned bootstrap current, that is carried by the passing (i.e., non-trapped) particles [1]. With a magnetic island present, the flattening of the pressure gradient creates a hole in the bootstrap current, and (for typical tokamak conditions) the resulting filamentation of current density leads to an amplification of the magnetic island. This is the neoclassical tearing mode instability [2-9]-a major concern for ITER because, if not controlled, it causes significant confinement degradation and can even terminate the plasma discharge in a disruption.
Neoclassical theory is well-developed to describe the physics of the trapped and passing particles when their trapped banana orbit widths, ρ bi and ρ be , are much smaller than the length scales of the system. The theory therefore provides a good description of the bootstrap current drive when the islands are much wider than the ion banana width. The NTM theory then predicts that all seed magnetic islands, however small, will grow to large amplitude with Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. a deleterious impact on confinement. However, experiments indicate that the plasma heals sufficiently small magnetic islands of half-width w ∼ few ρ bi [10]. This is precisely the regime where the conventional theory breaks down. Indeed, Ref. [11] demonstrated with particle-in-cell simulations that an ion density gradient is supported for islands in this regime, and the bootstrap current perturbation is suppressed. However, that work did not address the electron response. We find this introduces new physics and thus, to develop a quantitative understanding of the threshold phenomenon, it is necessary to understand (i) how the ions respond to small islands, (ii) the implications for the electron response via quasineutrality, and (iii) the consequences for the NTM drive. We address each of these in this Letter.
Our starting point is the drift kinetic model to describe the ion distribution function in a magnetized plasma, with electrostatic potential Φ. Assuming that the effect of the island on plasma parameters is localized to its vicinity, we work in the island rest frame and seek a steady state solution, neglecting any temporal fluctuation in fields (see Ref. [12] for the impact of turbulence). While finite ion Larmor radius effects could be included using a gyrokinetic approach, they are not essential for the physics we describe here. We therefore adopt the drift kinetic equation for the ion distribution function where v is the particle speed (k denoting a component parallel to the magnetic field, B), Þ is the combination of grad-B and curvature drifts, ω ci ¼ ZeB=m i , and Ze and m i are the ion charge and mass, respectively. On the right-hand side of Eq. (1) is the collision operator. Spatial derivatives are taken at constant kinetic energy where ψ is the poloidal flux and ψ s labels the rational surface where the safety factor, qðψ ¼ ψ s Þ¼m=n, with m and n being the poloidal and toroidal mode numbers of the magnetic island, respectively. Our other two spatial coordinates are the straight field line poloidal angle θ measuring the distance along an equilibrium magnetic field line, and helical angle ξ labeling the field lines at the rational surface. For a toroidally symmetric tokamak plasma, the canonical angular momentum, p ϕ ¼ðψ − ψ s Þ − Iv k =ω ci , is conserved during particle motion, where IðψÞ¼RB ϕ with R the major radius and B ϕ the toroidal component of the magnetic field. Exploiting this variable is key to reducing the dimensionality of the system. We solve Eq. (1) by writing the ion distribution Mi ð0Þþg i , where (0) indicates the quantity evaluated at the rational surface, the prime denotes a differential with respect to ψ, T i is the ion temperature, and is the ion thermal speed, n i is the ion density, and k B is Boltzmann's constant.
The magnetic field of a tokamak has a maximum, B max , on the inboard side. If λ < λ c ≡ B −1 max , the particles pass around the full extent of the flux surface, most of them deviating from it by a small amount ∼ϵρ θi , where ρ θi ¼ m i v thi =ZeB θ is the poloidal ion Larmor radius, and ϵ ¼ r=R is the inverse aspect ratio. If λ > λ c the particles are trapped to the region of weaker magnetic field, bouncing between the two points θ bAE along the field lines, where λBðθ ¼ θ bAE Þ¼1, and deviating from the flux surface by a larger amount: the ion banana width, ρ bi ∼ ϵ 1=2 ρ θi .I na tokamak, the system size is typically much greater than ρ θi , so we can introduce a small parameter, Δ ¼ ρ θi =r s , where r s is the minor radius of the rational surface where ψ ¼ ψ s . We consider the ion response to small islands with a width w ∼ ρ bi , and seek an asymptotic series solution to Eq. (1) by expanding in powers of Δ: Ordering w=r, ZeΦ=T i , and g ð0Þ i =F Mi ð0Þ all like Δ, the leading order contributions to Eq. (1) come from the freestreaming along the magnetic field lines, as well as the radial components of the grad-B and curvature drifts, which combine to give v k B assuming that collisions are OðΔÞ smaller than this freestreaming term. Integrating Eq. (2) shows that g i is independent of θ at fixed p ϕ , i.e., g ð0Þ i ðx; θ; ξ; λ;v; σÞ1 g ð0Þ i ðp ϕ ; ξ; λ;v; σÞ. This reduces the dimension of the problem, but the form ofḡ has the same form as that in Eq. (2). This term is eliminated by taking the average over θ at fixed p ϕ , ξ, λ, and v, which is equivalent to averaging over the particle orbits. For λ < λ c , we can integrate over a period in θ, imposing a periodic boundary condition to eliminate the term in g ð1Þ i . In the trapped region (λ > λ c ), we integrate between the bounce points θ bAE and sum over the two streams, σ ¼AE1, which then annihilates the term in g ð1Þ i due to continuity at each of the bounce points. Considering a large aspect ratio circular cross section tokamak, dropping terms of Oðϵ 2 Þ and smaller, and writing the magnetic field perturbation B 1 ¼ ∇ × ðA k bÞ, with RA k ¼ −w 2 q 0 =ð4qÞfðξÞ (q 0 ¼ dq=dψ), we arrive at our final equation forḡ where h…i θ denotes the orbit average described above, and F ðξÞ¼cos ξ. We have defined dimensionless variableŝ p ϕ ¼p ϕ =ψ s ,v jj ¼v jj =v thi , ν Ãi ¼ν i Rq=ðϵ 3=2 v thi Þ andŵ¼w=ψ s , together with the dimensionless drift frequencies Φ ¼ L q ðZe=T i ÞΦ, L −1 q ¼ð1=qÞdq=dx,ρ θi ¼ ρ θi =ψ s and Θ is the Heaviside function. Equation (3) is the solubility condition for g ð0Þ i . We solve Eq. (3) numerically for arbitrary values ofρ θi andŵ, employing the momentum-conserving model collision operator described in Ref. [13]. To determine the electrostatic potentialΦ we impose the quasineutrality condition, which requires a solution for the electron response. Because the electron orbit width is a factor ðm e =m i Þ 1=2 smaller than that of the ions, we adopt the small ρ θe =w solution described in Ref. [5]. To ensure that the collisions correctly account for momentum conservation, we use our numerical solutions for the ion flow in the electron collision operator. Figure 1 shows a color contour plot of the passing ion distribution function forŵ ¼ρ θi ¼ 0.02, L q =r s ¼ 1.0, λ=λ c ¼ 0.1,v ¼ 1.0, ν Ãi ¼ 0.01, and ϵ ¼ 0.1 (likewise for subsequent figures, unless otherwise stated). The island structure is clear in the color contours, but comparison with the magnetic island flux contours shows a shift in the contours of constant distribution function compared to the magnetic island. To understand this, consider the collisionless limit of Eq. (3) which, to leading order in Δ, can be reduced to the following form: This function S defines the stream lines for the ions. When the effect of Φ is ignored (which represents the effect of the E × B drift), one can show that the contours of constant S are identical to the magnetic island flux surfaces, but shifted by a fewρ θi . The result of Eq. (5) is that the ion distribution function now only depends on three variables in the low collision frequency limit: g ð0Þ i ðx; θ; ξ; λ;v; σÞ1 g ð0Þ i ðS; λ;v; σÞ. This is the toroidal generalization of the result from slab geometry [i.e., Eq. (8) of Ref. [14]]. The dependence on S, λ, and v can be derived by introducing collisions at next order to provide another constraint equation [15]. The contours of constant S are shown as the full curves in Fig 1, confirming that they coincide with the color contours of the distribution function. We refer to the constant S island structures as "drift islands." The shift of the drift island for σ ¼þ1 is equal and opposite to that for σ ¼ −1. In constructing the density, one sums over σ before integrating the distribution function over λ and v. Because the regions where the distribution function is flattened shift in opposite directions for σ ¼AE1, the distribution summed over σ supports substantial gradient inside the magnetic island whenŵ ∼ρ θi . This is a finite orbit width effect-not the well-known transport effect [4,[16][17][18]. For largeŵ ≫ρ θi , the shift in each direction is relatively small, and then the density is approximately flattened across the magnetic island, as expected (see Fig. 2). That finite orbit width effects support a pressure gradient in the ions is not new, and was also found numerically in Ref. [11]. Here we have shown that it is a consequence of the drift islands of passing particles (which will be captured in their numerical simulations also) and not, as suggested in Ref. [11], a result of trapped ions averaging over inside and outside the magnetic island as they intersect its separatrix. New physics also arises from the electron response, which was not retained in Ref. [11]. The strong parallel flows of the electrons tend to flatten their density across the magnetic island even for the small island width case. However, the electrostatic potential adjusts to ensure that the adiabatic part of the electron response provides the same gradients in the electron density as we see here for the ion density. This has profound consequences for NTM stability as we discuss shortly.
Let us first consider the ion flows. For them, the parallel streaming and E × B flows are expected to compete (e.g., see Refs. [5,13]), so it is particularly important that the electrostatic potential is derived consistently with the quasineutrality condition. From Fig. 3, it is clear for the smallρ θi =ŵ case that the pattern of the flow around the island follows the perturbed magnetic island geometry. On the other hand, the case forρ θi ¼ŵ is more complicated-there is a notable variation in the flow within a flux surface, with a broader peak in the profile a few island widths beyond the separatrix. This may provide an experimental test.
To calculate the impact on the NTM drive, we require the current perturbation. Combining our numerically derived ion flows with the analytic theory for the electron neoclassical flows [5], we calculate the current averaged over the magnetic island flux surfaces. The contribution to the island evolution is characterized by Δ 0 bs , given by [5] ψ s where J bs ¼ P j¼i;e Z j en j v thj hu kj i Ω and hÁÁ Ái Ω ¼ with Ω ¼ 2x 2 =ŵ 2 − cos ξ the perturbed flux function describing the island geometry. The flux surface integrals are taken at constant Ω.
The results for Δ 0 bs normalized to β θ as a function ofŵ are shown for a range ofρ θi in Fig. 4 where p is the plasma pressure). For largeŵ ≫ρ θi the result asymptotes to the value expected from previous analytic theories [2,3,5], represented by the dashed line. However, for small island widths we see that the impact of the shift of the drift islands is to reduce the bootstrap drive. The negative value of Δ 0 bs at the smallest island widths indicates that the effect of the current perturbation is to heal the island-a remarkable and unexpected result (not seen in the simulation results of Ref. [11], for example). For larger ρ θi , the peak value in Δ 0 bs decreases substantially, suppressing the bootstrap drive for the island growth. The critical island width, w c , where Δ 0 bs passes through zero, increases linearly with ρ θi : it can be fitted by w c ≃ 2.7ρ θi .
To understand the stabilization of small islands, we plot in Fig. 5 the individual ion and electron current density contributions to Δ 0 bs . Plotting ρ θi Δ 0 bs =β θ vs w=ρ θi , we find that all five ρ θi =r cases condense onto a universal set of curves for the ion and electron contributions. This is a consequence of the parallel flows being proportional to ρ θi;e , as predicted by analytic neoclassical theory. Notice that, as w → 0, the ion contribution to Δ 0 bs tends to zero, consistent with the bootstrap current being unperturbed in this limit (as found in Ref. [11]). The electron contribution to Δ 0 bs , on the other hand, is strongly negative at small w,so it is those that are responsible for healing small islands. As  the ion bootstrap current perturbation tends to zero with w, it is natural to assume the electron part will also, and so we postulate that the healing of the islands is a consequence of the electron response to the potential that is required to maintain quasineutrality.
In conclusion, we have presented a new drift kinetic theory for the response of ions to small magnetic islands in tokamak plasmas, and deduced some of the implications for the neoclassical tearing mode threshold physics. Neglecting cross field transport, we find that a consequence of the drifts is that the ion distribution function is not flattened across the magnetic island, but rather across a drift island that is shifted radially compared to the magnetic island. This shift is important for small islands comparable to the trapped ion banana width, in which case a pressure gradient is maintained inside the island, explaining previous simulation results [11]. This suppresses the bootstrap current drive for the NTM and the flows are then dominated by the electron physics, tending to heal a sufficiently small seed island. This new physics is important for a complete theory of the neoclassical tearing mode threshold and, in particular, for designing the NTM control system for ITER. Understanding the full implications of our theory for quantifying the NTM threshold will be the subject of future work, including an assessment of the accuracy of our analytic theory for the electron response employed here, the impact of finite ion Larmor radius, and the impact of finite island propagation frequency, including the ion polarization current physics. This work was supported by the UK Engineering and Physical Sciences Research Council, Grant No. EP/ N009363/1. Numerical calculations were performed using the ARCHER computing service through the Plasma HEC Consortium EPSRC Grant No. EP/L000237/1, as well as on EUROfusion High Performance Computer (Marconi-Fusion).