New insights into the generalized Rutherford equation for nonlinear neoclassical tearing mode growth from 2D reduced MHD simulations

Two dimensional reduced MHD simulations of neoclassical tearing mode growth and suppression by ECCD are performed. The perturbation of the bootstrap current density and the EC drive current density perturbation are assumed to be functions of the perturbed flux surfaces. In the case of ECCD, this implies that the applied power is flux surface averaged to obtain the EC driven current density distribution. The results are consistent with predictions from the generalized Rutherford equation using common expressions for Δbs′ ?> and ΔECCD′ ?>. These expressions are commonly perceived to describe only the effect on the tearing mode growth of the helical component of the respective current perturbation acting through the modification of Ohm’s law. Our results show that they describe in addition the effect of the poloidally averaged current density perturbation which acts through modification of the tearing mode stability index. Except for modulated ECCD, the largest contribution to the mode growth comes from this poloidally averaged current density perturbation.


Introduction
Although the generalized Rutherford equation for the non linear evolution of magnetic island widths is commonly used to model neoclassical tearing mode (NTM) growth [1,2], a numerical validation of its ability to predict the growth rate using the underlying magnetohydrodynamic (MHD) model is conspicuously lacking. The generalized Rutherford equa tion is obtained by averaging the current diffusion equation for the dominant Fourier harmonic of the helical magnetic flux perturbation over the island region, where resistivity is impor tant. This averaged solution is then matched to the linear, ideal MHD solution outside the magnetic island region. During the process of averaging, all information on the detailed structure of the island is discarded. In particular, possible asymmetries and contributions from higher Fourier harmonics of the magn etic flux perturbation are neglected. At the same time, the gen eralized Rutherford equation has been used successfully to describe the experimentally observed growth of NTMs and their stabilization by localized electron cyclotron current drive (ECCD). This has motivated the extensive use of the general ized Rutherford equation to establish requirements on electron cyclotron current drive systems for the suppression of NTMs in both current and future tokamaks like ITER. However, 2D Two dimensional reduced MHD simulations of neoclassical tearing mode growth and suppression by ECCD are performed. The perturbation of the bootstrap current density and the EC drive current density perturbation are assumed to be functions of the perturbed flux surfaces. In the case of ECCD, this implies that the applied power is flux surface averaged to obtain the EC driven current density distribution. The results are consistent with predictions from the generalized Rutherford equation using common expressions for bs ∆ ′ and ECCD ∆ ′ . These expressions are commonly perceived to describe only the effect on the tearing mode growth of the helical component of the respective current perturbation acting through the modification of Ohm's law. Our results show that they describe in addition the effect of the poloidally averaged current density perturbation which acts through modification of the tearing mode stability index. Except for modulated ECCD, the largest contribution to the mode growth comes from this poloidally averaged current density perturbation.
Keywords: neoclassical tearing modes (NTM), electron cyclotron current drive (ECCD), generalized Rutherford equation, 2D reduced MHD (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. and 3D (reduced) MHD simulations suggest that our under standing of the nonlinear tearing mode growth in terms of the Rutherford equation is incomplete [3,4]. Yu et al [3] observe that the stabilizing effect of an EC driven current with a depo sition width w cd smaller than the magnetic island width w, can be attributed in almost equal amounts to two effects. The first is due to the poloidally averaged EC driven current den sity. The second contribution comes from the helical pertur bation of the EC driven current density that results from the fast equilibration of the EC driven current density along the closed magnetic surfaces inside the magnetic island. In the parameter regime w w cd < , our current understanding of the generalized Rutherford equation predicts that stabilization should be almost entirely due to the helical component of the EC driven current density, while the effect of the poloidally averaged current density has been predicted to be exponen tially small [5].
In this paper, we present results of simulations of the sup pression of NTMs by ECCD with a 2D reduced MHD code [6] that shed new light on the generalized Rutherford equa tion. We highlight the relative roles of the poloidally averaged current density perturbation and the helical component of the current density perturbation. In the Rutherford theory these two terms act in two different ways: through modification of the equilibrium tearing mode stability index ∆′ and through modification of Ohm's law, respectively. We demonstrate that the commonly used expressions in the generalized Rutherford equation for the effect of helical current density perturbations [2,7,8] already include the effect of the corresponding poloi dally averaged current density perturbation. With the excep tion of the case where the ECCD is modulated in phase with the island rotation, the poloidally averaged current density perturbation is responsible for the larger contribution to the suppression of the magnetic island.

Analytical considerations
The generalized Rutherford equation for the nonlinear evol ution of the full width w of a magnetic island associated with a (neoclassical) tearing mode can be written in a generic form as [2,5,8] where g 1 is a geometric constant equal to g 1 = 0.82 [9], and η is the resistivity. The terms on the right hand side represent the tearing mode stability index ∆′ in case of the unperturbed equilibrium, the modification j i ( ) δ δ ∆′ of the tearing stability index as a consequence of the poloidal average of the cur rent density perturbations j i δ , and the effect of the helical current density perturbation inside the island region j i ∆ ′ δ . The latter two terms are summed over all possible current den sity perturbations including, in particular, the perturbation of the bootstrap current which results from the presence of the magnetic island (the main destabilizing term driving NTMs), and the EC driven current density applied for the suppression of NTMs.
When a current density perturbation is localized near the resonant surface of the tearing mode, the resulting modifica tion of the tearing mode stability index has been shown to be given by [10,11] where x r r s ≡ − is a radial coordinate measuring the distance from the resonant surface, ( ) eq is the leading order in the Taylor expansion of the equilibrium helical magnetic flux (where the ′ denotes the derivative with respect to x), and is the poloidal average of the cur rent perturbation (where ξ is the helical angle). The symbol P denotes that in case of a singular integrand at the resonant surface, the principal value of the integral is to be taken. This expression assumes that the island width is much smaller than the radial extension of the current density perturbation. In the case of finite island width, the integration interval should be taken over the exterior region up to the point where it is matched to the solution for the island interior. This matching is then commonly done In a tokamak with regular shear the second derivative of the equilibrium helical flux is negative when the positive toroidal angle is chosen in the direction of the plasma current. The general expression for the effect of the helical current perturbation in the magnetic island region is [2,7,8] is the part of the cur rent perturbation with the same helicity as the magnetic island. In equation (4) the borders of the integration domain assume that the matching to the exterior linear solution is performed outside a region covering both the island as well as the entire current perturbation.
Fast parallel transport along the closed magnetic field lines of a magnetic island results in a constant pressure in that region. As a consequence, the bootstrap current, which is driven by a radial pressure gradient, is annihilated inside the island. Substituting in equation (4) a constant j J bs bs δ = − inside the magnetic island and a zero perturbation outside, the effect of the helical perturbation of the bootstrap current is evaluated analytically to give [5] bs eq bs (5) Since the perturbation of the bootstrap current j bs δ is entirely localized inside the magnetic island, it is assumed that it does not perturb the stability index and any effect of the poloidally averaged current density on the mode growth is expected to be negligible. Note that this expression does not include the effects from the competition between parallel and perpendicular transport, that result in an only partial annihila tion of the bootstrap current inside small islands [12].
An EC driven current density is often represented by a Gaussian of amplitude J cd and full width w cd centered on the rel ative position x cd , i.e.
( ) Assuming perfect localization of the EC driven current (x 0 cd = ), the resulting modification of the tearing mode sta bility index is ECCD eq cd cd (6) In the presence of a finite size magnetic island, Bertelli et al [5] suggest that this change in the tearing mode stability index must be corrected by a multiplication with w w erfc / cd ( ) to account for the change in the integration domain (3). This assumes that the presence of the island does not significantly affect the poloidally averaged current density profile outside the island region. Inside the island, however, the EC driven current density distribution is profoundly changed: because of fast parallel transport along the magnetic surfaces, the cur rent density distribution will effectively be changed into a flux surface average of the original profile. We indicate this flux surface averaged current density by is the normalized flux surface label in the presence of the magnetic island. This flux surface averaged EC driven current density possesses a helical component, which according to equation (4) contributes an additional term in the Rutherford equation. We write this contrib ution in the form as given in [8] where F CD is a geometric function depending on the island size, the width and position of the EC driven current density profile, and the duty cycle D of a possible power modula tion. It is customary in the literature to evaluate this geometric function using the infinite integration domain of equation (4). When instead the interior and exterior solutions are matched at the edge of the island region, as in the case of finite size islands, the integration domain in equation (4) is limited to complementary to the integration domain (3) determining the contribution of the poloidally aver aged current density perturbation outside the island region to the tearing mode stability index. Figure 1 shows the results of a numerical evaluation of these two contributions from continu ously applied (CW) ECCD, which has a duty cycle D = 1, as obtained over these complementary integration domains. Also the sum of these two contributions is shown, which is seen to be almost indistinguishable from the value ∆ ′ ECCD of the contrib ution from only the helical current perturbation obtained by evaluating the radial integral over the infinite domain. For the latter the analytical fit to F CD as obtained in [8] is used: Similarly, using a partial integration and ξ The integrands of the radial integration of equations (9) and (10) are identical up to a contribution from the second Fourier harmonic of j i ⟨ ⟩( ) δ Ω , which is generally negligible. In conclusion, no matter how the two complementary radial integration domains are chosen, the modification of the tearing mode stability index and the effect from the helical current perturbation inside the island always sum approximately to ∆ ′ δj as evaluated over the infinite integra tion domain. Consequently, using this expression in the gen eralized Rutherford equation accounts for the full effect of a current density perturbation on the growth of a tearing mode, including both the modification of the tearing mode stability index due to the poloidally averaged current density as well as the effect of the helical component of the current density inside the island region.

Description of the 2D reduced MHD code
In 2D the dynamics of a tearing mode can be described by the reduced MHD equations for the helical magnetic flux ψ and the potential ϕ [9]: where η is the resistivity, ν the viscosity, j 2 ψ ≡ ∇ the current density, e Be B z zẑψ = ×∇ + the magnetic field, e v ẑ ϕ = ×∇ the velocity, and j ni δ represents all perturbations to noninduc tively driven currents (bootstrap current, ECCD, etc.). The constant electric field E is set to maintain the equilibrium cur rent density.
We have developed a numerical code that solves these equations with the specific purpose of validating the Rutherford equation. In line with the approximations under lying the Rutherford equation, we use an equilibrium helical magnetic flux given by the leading order of its Taylor expan sion around the resonant surface x = 0: In the same spirit, the radial boundary condition for each Fourier harmonic of the helical magnetic flux perturbation at x = −L and +L is determined by the step in its logarithmic deriva tive consistent with the value of its exterior stability index ∆ ′ k (given as an input to the code). These boundary conditions are enabled by the use of a mixed finite difference discretization in the radial direction x and Fourier decomposition in the peri odic angular direction ξ. Similarly, the boundary condition for the stream function of the dominant Fourier mode is obtained by matching to the linear ideal MHD solution in the exterior region. The code reproduces both the linear growth rate of an unstable tearing mode, as well as the Rutherford growth for finite size magnetic islands [6].
The effect of the bootstrap current perturbation is included in the code by a constant current density perturbation equal to J bs − inside the entire island region. This is then projected onto the set of Fourier harmonics in the calculation. A pos sible EC driven current density profile, which is assumed to be Gaussian in the radial direction, is averaged over the magnetic flux surfaces. The flux surface average is then projected onto the set of Fourier harmonics in the code. We have verified that the modification of the linear growth rate as reproduced by the code in the presence of a well localized EC driven current, is consistent with the modification of the tearing mode stability index as given in equation (6).

Results of numerical simulations
In order to study the effects from the bootstrap current per turbation and ECCD on the growth of magnetic islands, we perform simulations using the following parameters: in all runs, the equilibrium helical flux (normalized by 0 ρµ ) is 5  In order to set up an initial state for our simulations of NTM growth and suppression, we start a simulation from a helical magnetic flux perturbation which represents a 3 cm wide island. After a short phase in which the flow adjusts to the presence of the initial helical flux perturbation, the island shrinks in size at a constant rate consistent with the Rutherford equation. The simulation is stopped when an island size of about 0.5 cm is reached, which forms the ini tial state of subsequent simulations. We simulate a NTM that is seeded with an initial width of the magnetic island of w 0.5 cm = , and is subsequently suppressed by ECCD (see figure 2). The perturbation to the bootstrap current inside the island is j 6630 bs δ = − . For this value of the bootstrap current, the Rutherford equation predicts a saturated island size of approximately 9 cm. When the NTM reaches an island size of 3 cm, however, ECCD is switched on with a maximum driven current density of J 15 000 cd = . The ECCD is centered exactly at the resonance surface x 0 cd = with a Gaussian profile width of w 1 cd = cm. Two cases are simulated: one for CW ECCD and the other for modulated ECCD with a duty cycle of 50% centered around the Opoint phase of the magnetic island. The maximum driven current density is identical for the CW and modulated cases, resulting in a factor of two reduction of the total driven current in the case of modulated ECCD. Because of the way we have implemented the bootstrap cur rent perturbation, it is not reduced for small island sizes. Thus full suppression of the mode by ECCD is not expected in our simulations. The simulation is stopped when the island reaches a size of about 0.5 cm.
From the initial growth phase of the NTM, we obtain an estimate of the contribution bs ∆ ′ from the bootstrap current perturbation to the tearing mode growth. We compare this with the analytical prediction for this contribution as given in equation (5). In order to study the contributions of individual harmonics of the bootstrap current perturbation, the simula tion is repeated taking into account either only the poloidally averaged current density j k bs, 0 δ = or only its helical component j k bs, 1 δ = . The results are shown in figure 3. The code results are in excellent agreement with the analytical prediction. However, we must note here, that while the analytical predic tion appears to be based solely on the helical component of the current perturbation, the poloidal average rather than the helical perturbation provides the largest contribution to bs ∆ ′ in the simulations. Because the k = 0 poloidally averaged per turbation acts through a modification of the stability index, it requires a finite current diffusion time to take effect. This explains the deviation from the analytical curve at small island sizes, which represents the time frame immediately after the seeding of the mode. In contrast, the effect of the k = 1 helical perturbation is instantaneous.
During the phase of island suppression by ECCD, the island size samples a large range of the function F CD covering both island sizes larger and smaller than the ECCD profile width. Assuming that the full effect of the ECCD is represented by ∆ ′ ECCD as given in equation (7), an estimate of F CD is obtained from the rate of change of the island width observed in the simulations. In figure 4 the results for (a) CW ECCD and (b) modulated ECCD are compared with the analytical expecta tion obtained from the fit functions as given by De Lazzari et al [8]. Once more the agreement between the code results and the analytical prediction is excellent. Also results obtained by keeping only the k = 0 or k = 1 contribution from j ECCD ⟨ ⟩ in the code are shown. As for the bootstrap current, the effect of CW ECCD must be ascribed mostly to the change in the equi librium stability generated by the poloidally averaged cur rent density perturbation, even for the cases where the island size is significantly larger than the ECCD profile width. In fact, the results for k = 1 require the use of ten times the CW EC driven current density to obtain a sufficiently stabilizing effect. The modulated ECCD with a duty cycle of 50% is seen to be more efficient then CW ECCD over the entire range of parameters covered by the simulations. The advantage of modulated ECCD over unmodulated ECCD is particularly large at small island sizes consistent with previous analysis [13] and 3D MHD numerical simulations [14]. In this case the larger contribution comes from the k = 1 helical component of the EC driven current density.

Discussion and conclusions
The results of the simulations clearly demonstrate that the gen eralized Rutherford equation provides an excellent description of the nonlinear growth of NTMs. The effect of the bootstrap current perturbation and the ECCD is shown to be fully con tained within the expressions for bs ∆ ′ (5) [5] and ECCD ∆ ′ (7) [8]. However, contrary to the common assumption in the literature that these expressions contain solely the effect of the helical current perturbation, we have shown that they contain both the effect of the helical part of the current density perturbation as well as the modification of the equilibrium tearing mode sta bility index from the poloidally averaged current perturbation. With the exception of the case of modulated ECCD, we have demonstrated numerically that the effect of the poloidally averaged current density perturbation is dominant.
In the theory giving rise to the Rutherford equation, the terms representing the contributions to the growth of the mode arising from the poloidally averaged k = 0 current den sity perturbation and its k = 1 helical component come from   (7), respectively, as evaluated over the infinite integration domain [5,8]. We conclude that these expres sions thus include both the modification the tearing mode stability index and the effect of the helical current density perturbation.
The two effects can be distinguished in our 2D reduced MHD simulations. The effect of the poloidally averaged cur rent density perturbation on the modified stability index is only realized after a current diffusion time over the width of the driven current density profile. In contrast, the effect of the helical current perturbation acting through the modi fication of Ohm's law inside the magnetic island is instanta neous. This is clearly seen in the results of the 2D reduced MHD code. Including an additional term ( ) δ∆′ j ECCD in the generalized Rutherford equation, as is sometimes done [2,5], would result in double counting of the effect of the poloidally averaged EC driven current density on the tearing mode stability index.
In order to compare the relative contributions from the poloidally averaged k = 0 component and the helical k = 1 component to analytical expectations, the simulation results for CW ECCD in figure 4(a) can be compared to the ana lytical results for these contributions from CW ECCD based on a matching of the exterior and interior solutions at the edge of the magnetic island, x w/2 = ± , as presented in figure 1. The analytical result clearly underestimates the contribution from the modification of the stability caused the poloidally averaged k = 0 component of the current density perturba tion. In the case of the bootstrap current perturbation which is completely localized inside the island the usual matching at x w/2 = ± , would predict a negligible contribution through the modification of the stability index from the k = 0 comp onent. This is again in clear contrast to the simulations. In conclusion, all simulation results indicate a matching of exte rior and interior solutions at a radial position well inside the island region. The ECCD suppression rate as obtained from 2D reduced MHD simulations (full curves). Two cases are shown: (a) CW ECCD and (b) modulated ECCD with a duty cylcle of 50%. The results are presented in terms of an equivalent function F CD assuming that the full effect of the ECCD is contained in equation (7). The dotted curves give the analytical expectation for F CD as obtained from the fit functions given in [8]: for CW ECCD