Effect of Nonlinear Energy Transport on Neoclassical Tearing Mode Stability in Tokamak Plasmas

An investigation is made into the effect of the reduction in anomalous perpendicular electron heat transport inside the separatrix of a magnetic island chain associated with a neoclassical tearing mode in a tokamak plasma, due to the flattening of the electron temperature profile in this region, on the overall stability of the mode. The onset of the neoclassical tearing mode is governed by the ratio of the divergences of the parallel and perpendicular electron heat fluxes in the vicinity of the island chain. By increasing the degree of transport reduction, the onset of the mode, as the divergence ratio is gradually increased, can be made more and more abrupt. Eventually, when the degree of transport reduction passes a certain critical value, the onset of the neoclassical tearing mode becomes discontinuous. In other words, when some critical value of the divergence ratio is reached, there is a sudden bifurcation to a branch of neoclassical tearing mode solutions. Moreover, once this bifurcation has been triggered, the divergence ratio must reduced by a substantial factor to trigger the inverse bifurcation.


I. INTRODUCTION
Neoclassical tearing modes are large-scale magnetohydrodynamical instabilities that cause the axisymmetric, toroidally-nested, magnetic flux-surfaces of a tokamak plasma to reconnect to form helical magnetic island structures on low mode-number rational magnetic flux surfaces. 1 Island formation leads to a degradation of plasma energy confinement. 2 Indeed, the confinement degradation associated with neoclassical tearing modes constitutes a major impediment to the development of effective operating scenarios in present-day and future tokamak experiments. 3 Neoclassical tearing modes are driven by the flattening of the temperature and density profiles within the magnetic separatrix of the associated island chain, leading to the suppression of the neoclassical bootstrap current in this region, which has a destabilizing effect on the mode. 4 The degree of flattening of a given profile (i.e., either the density, electron temperature, or ion temperature profile) within the island separatrix depends on the ratio of the associated perpendicular (to the magnetic field) and parallel transport coefficients. 5 The dominant contribution to the perpendicular transport in tokamak plasmas comes from small-scale drift-wave turbulence, driven by plasma density and temperature gradients. 1 The fact that the density and temperature profiles are flattened within the magnetic separatrix of a magnetic island chain implies a substantial reduction in the associated perpendicular transport coefficients in this region. Such a reduction has been observed in gyrokinetic simulations, [6][7][8][9] as well as in experiments. [10][11][12][13][14] A strong reduction in perpendicular transport within the magnetic separatrix calls into question the conventional analytic theory of neoclassical tearing modes in which the perpendicular transport coefficients are assumed to spatially uniform in the island region. 5 The aim of this paper is to investigate the effect of the reduction in perpendicular transport inside the separatrix of a neoclassical magnetic island chain, due to profile flattening in this region, on the overall stability of the mode. For the sake of simplicity, we shall only consider the influence of the flattening of the electron temperature profile on mode stability. However, the analysis contained in this paper could be generalized, in a fairly straightforward manner, to take into account the influence of the flattening of the ion temperature and density profiles.

A. Fundamental Definitions
Consider a large aspect-ratio, low-β, circular cross-section, tokamak plasma equilibrium.
Let us adopt a right-handed cylindrical coordinate system (r, θ, z) whose symmetry axis (r = 0) coincides with the magnetic axis of the plasma. The system is assumed to be periodic in the z-direction with period 2π R 0 , where R 0 is the simulated major plasma radius. It is helpful to define the simulated toroidal angle ϕ = z/R 0 . The coordinate r serves as a label for the unperturbed (by the tearing mode) magnetic flux-surfaces. Let the equilibrium toroidal magnetic field, B z , and the equilibrium toroidal plasma current both run in the +z direction.
Suppose that a neoclassical tearing mode generates a helical magnetic island chain, with m θ poloidal periods, and n ϕ toroidal periods, that is embedded in the aforementioned plasma.
The island chain is assumed to be radially localized in the vicinity of its associated rational surface, minor radius r s , which is defined as the unperturbed magnetic flux-surface at which q(r s ) = m θ /n ϕ . Here, q(r) is the safety-factor profile (which is assumed to be a monotonically increasing function of r). Let the full radial width of the island chain's magnetic separatrix be W . In the following, it is assumed that ǫ s ≡ r s /R 0 ≪ 1 and W/r s ≪ 1.
It is convenient to employ a frame of reference that co-rotates with the magnetic island chain. All fields are assumed to depend (spatially) only on the radial coordinate r and the helical angle ζ = m θ θ − n ϕ ϕ. Let k θ = m θ /r s , q s = m θ /n ϕ , and s s = d ln q/d ln r| rs .
The magnetic shear length at the rational surface is defined L s = R 0 q s /s s . Moreover, the unperturbed (by the magnetic island) electron temperature gradient scale-length at the rational surface takes the form L T = −1/(d ln T 0 /dr) rs , where T 0 (r) is the unperturbed electron temperature profile. In the following, it is assumed that L T > 0, as is generally the case in conventional tokamak plasmas. 1 The helical magnetic flux is defined where x = r − r s , and the magnetic field perturbation associated with the tearing mode is written δB = ∇ × (δχ e z ). It is easily demonstrated that B · ∇χ = 0, where B is the total magnetic field. 15 Hence, χ is a magnetic flux-surface label. It is helpful to introduce the normalized helical magnetic flux, ψ = (L s /B z w 2 ) χ, where w = W/4. The normalized flux in the vicinity of the rational surface is assumed to take the form 15 where X = x/w. As is well-known, the contours of ψ map out a symmetric (with respect to X = 0), constant-ψ, 16 magnetic island chain whose O-points lie at ζ = π, X = 0, and ψ = −1, and whose X-points lie at ζ = 0, 2π, X = 0, and ψ = +1. The chain's magnetic separatrix corresponds to ψ = +1, the region inside the separatrix to −1 ≤ ψ < 1, and the region outside the separatrix to ψ > 1. The full radial width of the separatrix (in X) is 4.
Finally, the electron temperature profile in the vicinity of the rational surface is written where T s = T 0 (r s ), and Note that δT (X, ζ) is an odd function of X. In the following, it is assumed that w/L T ≪ 1.

B. Electron Energy Conservation Equation
The steady-state electron temperature profile in the vicinity of the island chain is governed by the following well-known electron energy conservation equation: 5,17 where and Here, κ ⊥ and κ are the perpendicular (to the magnetic field) and parallel electron thermal conductivities, respectively. The first term on the right-hand side of Eq. and ∂ 2 δT /∂X 2 are both O(1) in our normalization scheme, the ratio of the divergences of the parallel and perpendicular heat fluxes is effectively measured by (W/W c ) 4 .] The quantity W c is the critical island width above which the former term dominates the latter, causing the temperature profile to flatten within the island separatrix. 5 In writing Eq. (5), we have neglected any localized sources or sinks of heat in the island region. We have also assumed that κ ⊥ and κ are spatially uniform in the vicinity of the rational surface. The latter assumption is relaxed in Sect. III

C. Narrow-Island Limit
Consider the so-called narrow-island limit in which W ≪ W c . 5 Let Equation (5) transforms to give We can write where subject to the boundary conditions T 1 (0, ζ) = 0, and T 1 → 0 as |Y | → ∞. Note that the solution (10) automatically satisfies the boundary condition (4). It follows that where f (p) is the well-behaved solution of that satisfies f (0) = 0, and f → 0 as |p| → ∞. Note that f (−p) = −f (p). Hence, in the narrow-island limit, 5 D. Wide-Island Limit Consider the so-called wide-island limit in which W ≫ W c . 5 We can write whereT andT are both O(1), and T = 0.
Equations (5) and (15) can be combined to give The flux-surface average of the previous equation yields which implies that d dψ The previous equation can be integrated to givē which satisfies the boundary condition (4). Hence, in the wide-island limit, 5 (24)

E. Modified Rutherford Equation
The temporal evolution of the island width is governed by the so-called modified Rutherford equation, which takes the form 4,5,15 where Here, τ R = µ 0 r 2 s /η(r s ) is the resistive evolution timescale at the rational surface, and η(r) is the unperturbed plasma resistivity profile. Moreover, ∆ ′ < 0 is the standard linear tearing stability index. 16 z , and n s is the unperturbed electron number density at the rational surface. The second term on the right-hand side of Eq. (25) parameterizes the destabilizing influence of the perturbed bootstrap current. 4,5 Note that, in this paper, for the sake of simplicity, we have employed the so-called lowest-order asymptotic matching scheme described in Ref. 18. This accounts for the absence of higher-order island saturation terms in Eq. (25).

A. Introduction
In conventional tokamak plasmas, the dominant contribution to the perpendicular electron thermal conductivity, κ ⊥ , comes from small-scale drift-wave turbulence driven by electron temperature gradients. 1 The fact that the electron temperature gradient is flattened within the magnetic separatrix of a sufficiently wide magnetic island chain implies a substantial reduction in the perpendicular electron thermal conductivity in this region. There is clear experimental evidence that this is indeed the case. 11,12,14 In particular, Ref. 14 reports a reduction in κ ⊥ at the O-point of the magnetic island chain associated with a typical neoclassical tearing mode by 1 to 2 orders of magnitude. Obviously, such a strong reduction in κ ⊥ within the magnetic separatrix calls into question the conventional analytic model of neoclassical tearing modes, described in Sect. II, in which κ ⊥ is assumed to spatially uniform in the vicinity of the rational surface.

B. Nonuniform Perpendicular Electron Conductivity Model
As a first attempt to model the reduction in κ ⊥ due to temperature flattening within the magnetic separatrix of a neoclassical island chain, let us write where κ ⊥ 1 and κ ⊥ 0 are spatial constants, with κ ⊥ 1 ≤ κ ⊥ 0 . Since the mean temperature gradient outside the separatrix of a neoclassical magnetic island chain is similar in magnitude to the equilibrium temperature gradient [see Eq. (42)], it is reasonable to assume that κ ⊥ 0 is equal to the local (to the rational surface) perpendicular electron thermal conductivity in the absence of an island chain. Let be the critical island widths outside and inside the separatrix, respectively. Likewise, let measure the ratios of the divergences of the parallel and perpendicular electron heat fluxes outside and inside the separatrix, respectively. Finally, let the parameter measure the relative reduction of perpendicular electron heat transport within the island separatrix.
Let us adopt the following simple model: where 0 < δ ≤ 1. According to this model, the degree of perpendicular transport reduction within the separatrix is controlled by the parameter ξ 1 , which measures ratio of the divergences of the parallel and perpendicular electron heat fluxes inside the separatrix. (See Sects. II C and II D.) If ξ 1 is much less than unity then there is no temperature flattening within the separatrix, which implies that λ = 1 (i.e., there is no reduction in transport). On the other hand, if ξ 1 is much greater than unity then the temperature profile is completely flattened inside the separatrix, and the transport is reduced by some factor δ (say). The previous formula is designed to interpolate smoothly between these two extremes as ξ 1 varies.
Equations (33) and (34) can be combined to give If follows that δ ≤ λ ≤ 1, with ξ 0 = 0 when λ = 1, and ξ 0 → ∞ as λ → δ. It is easily demonstrated that the function ξ 0 (λ) has a point of inflection when δ = δ crit = 1/(1 + e 2 ) = 0.1192. This point corresponds to ξ 0 = 4 δ crit = 0.4768 and λ = 2 δ crit = 0.2384. Figure 1 shows the perpendicular electron transport reduction parameter, λ, plotted as a function of the ratio of the divergences of the parallel to perpendicular electron heat fluxes outside the island separatrix, ξ 0 , for various values of the maximum transport reduction parameter, δ. It can be seen that if δ > δ crit then the ξ 0 -λ curves are such that dξ 0 /dλ < 0 for δ ≤ λ ≤ 1. This implies that λ decreases smoothly and continuously as ξ 0 increases, and vice versa. We shall refer to these solutions as continuous solutions of Eq. is a bifurcation to the small-temperature-gradient solution branch. We shall refer to this bifurcation as the temperature-gradient-flattening bifurcation, because it is characterized by a sudden decrease in the transport ratio parameter, λ, which implies a sudden decrease in the electron temperature gradient within the island separatrix. Once on the small-temperaturegradient solution branch, the control parameter ξ 0 must be reduced significantly in order to trigger a bifurcation back to the large-temperature-gradient solution branch. We shall refer to this bifurcation as the temperature-gradient-restoring bifuration, because it is characterized by a sudden increase in the transport ratio parameter, λ, which implies a sudden increase in the electron temperature gradient within the island separatrix Figure 3 shows the critical values of the control parameter ξ 0 below and above which a temperature-gradient-flattening and a temperature-gradient-restoring bifurcation, respectively, are triggered, plotted as a function of δ/δ crit . Figure 4 shows the extents of the various solution branches (i.e., the continuous, largetemperature-gradient, small-temperature-gradient, and inaccessible branches) plotted in ξ 0ξ 1 space. It is clear that the large-temperature-gradient solution branch is characterized by ξ 0 ≪ 1 and ξ 1 < ∼ 1. In other words, the region outside the island separatrix lies in the narrowisland limit, W ≪ W c 0 , whereas that inside the separatrix lies in the narrow/intermediate (31) and (32).] This implies weak to moderate flattening of the temperature gradient within the separatrix. On the other hand, the small-temperarturegradient solution branch is characterized by ξ 0 ≪ 1 and ξ 1 ≫ 1. In other words, the region outside the island separatrix lies in the narrow-island limit, W ≪ W c 0 , whereas that inside the separatrix lies in the wide-island limit, W ≫ W c 1 . This implies strong flattening of the temperature gradient within the separatrix. Figure 4 suggests that bifurcated solutions of Eq. (35) occur because it is possible for the regions inside and outside the island separatrix to lie in opposite asymptotic limits (the two possible limits being the wide-island and the narrow-island limits). Obviously, this is not possible in the conventional model in which κ ⊥ is taken to be spatially uniform in the island region.
Finally, according to our simple model, the critical value of the maximum transport reduction parameter, δ, below which bifurcated solutions of the electron energy transport equation occur is 0.1192. As we have seen, there is experimental evidence for a transport reduction within the separatrix of a neoclassical island chain by between 1 and 2 orders of magnitude. 14 According to our model, such a reduction would be large enough to generate bifurcated solutions.

D. Evaluation of Integrals
According to Eqs. where Here, use has been made of the easily proved result Let ψ = 2 k 2 − 1. It follows that dψ = 4 k dk. In the region 0 ≤ k ≤ 1, we can write On the other hand, in the region k > 1, we can write Here, it is assumed that A is an even function of X. Let It follows from Eqs. (50)-(57) that in the region 0 ≤ k ≤ 1, On the other hand, in the region k > 1, Here, are complete elliptic integrals. 19 Hence, according to Eqs. (45)-(48) and (58)-(60), Thus, Eqs. (43) and (44) yield respectively.

E. Destabilizing Effect of Perturbed Bootstrap Current
The dimensionless parameter G 2 , appearing in the modified Rutherford equation, (25), measures the destabilizing influence of the perturbed bootstrap current. Figure 5 shows G 2 plotted as a function of the so-called neoclassical tearing mode control parameter, increases from a value much less than unity, we start off on the large-temperature-gradient solution branch, and the bootstrap destabilization parameter, G 2 , increases smoothly and monotonically from a small value. However, when a critical value of ξ 0 is reached, there is a gradient-flattening-bifurcation to the small-temperature-gradient solution branch. This bifurcation is accompanied by a sudden increase in G 2 to a value close to its asymptotic limit 6.140. Once on the small-temperature-gradient solution branch, ξ 0 must be decreased by a significant amount before a gradient-restoring-bifurcation to the large-temperature-gradient solution branch is triggered. Moreover, the gradient-restoring-bifurcation is accompanied by a very large reduction in G 2 .

F. Long Mean-Free-Path Effects
The parallel electron thermal conductivity takes the form 17 in a collisional plasma, where n e is the electron number density, v e is the electron themal velocity, and λ e is the electron mean-free-path. However, in a conventional tokamak plasma the mean-free-path λ e typically exceeds the parallel (to the magnetic field) wavelength λ of low-mode-number helical perturbations. Under these circumstances, the simple-minded application of Eq. (75) yields unphysically large parallel heat fluxes. The parallel conductivity in the physically-relevant long-mean-free-path limit (λ e ≪ λ ) can be crudely estimated as 5,20 κ ∼ n e v e λ , which is equivalent to replacing parallel conduction by parallel convection in the electron energy conservation equation, (5). For a magnetic island of full radial width W , the typical value of λ is n ϕ s s w/R 0 . Hence, in the long-mean-free-path limit, the expression for the neoclassical tearing mode control parameter (74) is replaced by where κ ′ = n ϕ n e v e s s /R 0 , and n e and v e are evaluated at the rational surface.

IV. SUMMARY AND DISCUSSION
In this paper, we have investigated the effect of the reduction in anomalous perpendicular electron heat transport inside the separatrix of a magnetic island chain associated with a neoclassical tearing mode in a tokamak plasma, due to the flattening of the electron temperature profile in this region, on the overall stability of the mode. Our model (which is described in Sect. III) is fairly crude, in that the perpendicular electron thermal conductivity, is characterized by relatively weak flattening of the electron temperature profile within the island separatrix, and a consequent relatively small value (i.e., significantly smaller than the asymptotic limit 6.140) of the bootstrap destabilization parameter, G 2 . On the other hand, the small-temperature-gradient branch is characterized by almost complete flattening of the electron temperature profile within the island separatrix. Consequently, the bootstrap destabilization parameter, G 2 , takes a value close to the asymptotic limit 6.140 on this solution branch. The large-temperature-gradient and small-temperature-gradient solution branches are separated by a dynamically inaccessible branch of solutions. As the control parameter, ξ 0 , increases from a value much less than unity, the system starts off on the large-temperature-gradient solution branch, and the bootstrap destabilization parameter, G 2 , increases smoothly and monotonically from a small value. However, when a critical value of ξ 0 is reached, there is a gradient-flattening-bifurcation to the small-temperaturegradient solution branch. (See Fig. 6.) This bifurcation is accompanied by a sudden increase in G 2 to a value close to its asymptotic limit 6.140. Once on the small-temperature-gradient solution branch, ξ 0 must be decreased by a significant amount before a gradient-restoringbifurcation to the large-temperature-gradient solution branch is triggered. Moreover, the gradient-restoring-bifurcation is accompanied by a very large reduction in G 2 . (See Section III E.) The behavior described in the preceding paragraph points to the disturbing possibility that a neoclassical tearing mode in a tokamak plasma could become essentially selfsustaining. In other words, once the mode is triggered, and the electron temperature profile is flattened within the island separatrix, the consequent substantial reduction in the perpendicular thermal conductivity in this region reinforces the temperature flattening, making it very difficult to remove the mode from the plasma.