The VMEC nonlinear ideal MHD equilibrium code [S. P. Hirshman and J. C. Whitson, Phys. Fluids **26**, 3553 (1983)] is compared against analytic linear ideal MHD theory in a screw-pinch-like configuration. The focus of such analysis is to verify the ideal MHD response at magnetic surfaces which possess magnetic transform (*ι*) which is resonant with spectral values of the perturbed boundary harmonics. A large aspect ratio circular cross section zero-beta equilibrium is considered. This equilibrium possess a rational surface with safety factor *q* = 2 at a normalized flux value of 0.5. A small resonant boundary perturbation is introduced, exciting a response at the resonant rational surface. The code is found to capture the plasma response as predicted by a newly developed analytic theory that ensures the existence of nested flux surfaces by allowing for a jump in rotational transform ($\iota =1/q$). The VMEC code satisfactorily reproduces these theoretical results without the necessity of an explicit transform discontinuity ($\Delta \iota $) at the rational surface. It is found that the response across the rational surfaces depends upon both radial grid resolution and local shear ($d\iota /d\Phi $, where *ι* is the rotational transform and $\Phi $ the enclosed toroidal flux). Calculations of an implicit $\Delta \iota $ suggest that it does not arise due to numerical artifacts (attributed to radial finite differences in VMEC) or existence conditions for flux surfaces as predicted by linear theory (minimum values of $\Delta \iota $). Scans of the rotational transform profile indicate that for experimentally relevant levels of transform shear the response becomes increasing localised. Careful examination of a large experimental tokamak equilibrium, with applied resonant fields, indicates that this shielding response is present, suggesting the phenomena is not limited to this verification exercise.

## I. INTRODUCTION

A screw-pinch with resonant boundary perturbation provides a means to validate the plasma response model of the VMEC 3D equilibrium code.^{1} While such a simple problem may lack experimental relevance, it allows us to clearly examine the equilibrium response to an applied resonant boundary perturbation at a magnetic surface with rational transform (those whose rotational transform can be expressed as *p*/*q*, where *p* and *q* are integers). Previous work by Newcomb suggests that a non-axisymmetric delta current will form on the rational surface to completely shield out the perturbed harmonic in the region interior to the surface.^{2} However, this results in overlapping of flux surfaces, a solution VMEC explicitly excludes through its magnetic field representation. Thus, a direct comparison between VMEC and the classical Newcomb's solution is ill-posed. Recent work by Loizu *et al.* indicates that solutions to Newcomb's equation, which do not cause flux surfaces to overlap, also exist.^{3} Comparisons between VMEC and these linear solutions are possible (avoiding the ill-posed comparison). In particular, we can better gauge how accurately VMEC is resolving the plasma response at resonant rational surfaces. Such analysis should also help us better understand discrepancies between benchmarking efforts focused on calculating the DIII-D tokamak experimental equilibria where resonant boundary perturbations are applied.

Recent efforts to benchmark various 3D equilibrium codes on the DIII-D tokamak have uncovered discrepancies between linear and non-linear ideal MHD codes. Initial work focused on examining the plasma response to applied *n* = 3 fields in DIII-D.^{4} In particular, the plasma response near magnetic surfaces with low order rotational transform was investigated. It was found that the VMEC non-linear 3D ideal MHD equilibrium code differed significantly from various linear codes. This was despite a plasma response, near rational surfaces, that scaled linearly with the plasma boundary deformation (validating the linear limit of the problem). However, in this work, it was noted that the nested flux surface constraint was violated (over the pedestal region) in the linear codes at relatively modest levels of external perturbation.

A second benchmarking effort on the DIII-D device examined the plasma response as measured at the vessel wall.^{5} In this work, the newly upgraded 3D magnetic diagnostic arrays were used to compare the measured plasma response against simulated results. Here, an *n* = 1 perturbation was applied using the DIII-D in-vessel coil set. The results suggested that while the linear codes and non-linear codes differed in their response at rational surfaces they agreed regarding the measured diagnostic response at the plasma wall. This served as a validation that VMEC could be used for 3D equilibrium reconstruction.^{6–9} It was suggested that a simplified model be examined to better understand the source of disagreement between the linear and non-linear ideal plasma response at a rational surface.

The simplest model to consider was a screw-pinch. Calculations performed by linear codes indicated the presence of a delta current at a resonant rational surface ($q=1/\iota =2$). Such a result is consistent with the classical Newcomb's solution to this problem.^{2} This current completely shields out the perturbation inside the rational surface, leading to a step-like response in the displacement of the flux surfaces. However, closer examination showed that arbitrarily close to the rational surface, flux surfaces were overlapping. The nested flux surface constraint could only be preserved for such models in the limit of zero external perturbation. The non-linear VMEC code, on the other hand, explicitly enforces the nested flux surface constraint. This demonstrates that, while Newcomb's solution could be used as a verification exercise for linear codes, comparison with codes like VMEC was dubious.

Work to reassess linear theory, focusing on the inclusion of the nested flux surface constraint, has recently been conducted by Loizu *et al*.^{3} It was determined that, in order to satisfy the nested flux surface constraint across a resonant rational surface, the initially unperturbed equilibrium must possess an axisymmetric sheet current at the resonant surface. As a consequence, the rotational transform jumps across the rational value (by $\Delta \iota $) and the surface no longer leads to a singularity. This sheet current was found to have a minimum value necessary to preserve the nested flux surface constraint. It was also found that despite the presence of a non-axisymmetric *δ*-current, on the resonant surface, the perturbation penetrates inside the resonant surface and a non-axisymmetric current is also established within the plasma volume. Moreover, for this class of MHD equilibria with discontinuous transform, an exact verification with a non-linear ideal MHD model was carried out (using the SPEC code^{10}). This verification exercise showed excellent convergence with respect to linear theory. It is now reasonable to ask how the VMEC code compares to these results.

This paper examines the VMEC screw-pinch-like equilibrium from the perspective of this newly established analytic model. Section II describes the configuration studied in VMEC and parameter scans that were performed. In Sec. III, the plasma response across the resonant surface is examined. Here, the possibility of an implicit axisymmetric delta current in VMEC is examined. A summary of the results is presented in the discussion with focus on the implication for more realistic geometries. The non-axisymmetric current density of a DIII-D equilibrium is also considered to highlight that these response features are not unique to the screw-pinch problem.

## II. METHOD

A perturbed large aspect ratio circular cross-section equilibrium serves as a basis for verification of the VMEC plasma response against analytic theory. There is an expectation for solutions to differ from that of Newcomb^{2} and linear codes because VMEC imposes the nested flux surface constraint everywhere. The linear codes have a step function discontinuity in their displacements at rational surfaces, which translates into an overlapping of surfaces.

The boundary perturbation and plasma profiles were chosen carefully so as to focus on the response at a *q* = 2 surface. The plasma pressure in these calculations was taken to be identically zero. This was done in order to suppress the $1/x$ type singularities associated with pressure gradients across rational surfaces.^{11,12} These singularities are associated with the resonant harmonic of the parallel current near the rational surface. An *ι* profile of the form $\iota (s)=0.6\u22120.2s$, where $s=\Phi /\Phi edge$ is the normalised toroidal flux, was chosen ($\Phi edge=6.28$ Wb). Figure 1 depicts the equilibrium configuration and *ι* profile. Such a profile attempts to avoid other resonant low-order rationals and interactions between multiple rational surfaces. In cylindrical geometry, the poloidal and toroidal Fourier modes do not couple. However, care must be taken with VMEC, as the code can only approximate cylindrical geometry (large aspect ratio, large field periodicity).

A boundary perturbation in minor radius (*ρ*) of the form $\delta \rho (n=1,m=2)$ was chosen so as to be resonant with the *q* = 2 surface. Here, *n* refers to the toroidal mode numbers and *m* the poloidal mode number. The VMEC code uses toroidal coordinates to define flux surfaces such that toroidal R and Z coordinates tracing out a flux surface have the form

where *s* is a radial flux surface label, *θ* is the poloidal angle, *ζ* is the toroidal angle, and *N _{FP}* is the fundamental periodicity of the problem. The values

*N*and

*M*define the truncation of the mode spectrum.

A transformation from the *ρ* minor radius representation to the toroidal coordinate representation must be performed. The cylindrical *ρ* representation of a surface can be written as

where *a* is the minor radius of the surface, $\delta \rho $ is the perturbed minor radius, and the subscripts on mode numbers are used to delineate toroidal from minor radius indices. Substituting into the relations $R=R00+\rho \u2009cos\u2009\theta $ and $Z=\rho \u2009sin\u2009\theta $, we find

where *R*_{00} is the major radius of the toroidal system. From these equations, the toroidal harmonics $R01=Z01=a,\u2009R1,3=R1,1=\delta \rho 1,2/2$ and $Z1,3=\u2212Z1,1=\delta \rho 1,2/2$ can be determined. More generally, for a given perturbed cylindrical harmonic $\delta \rho (n1,m1)$, the equivalent perturbed toroidal harmonics are $Rn1,m1+1=Rn1,m1\u22121=\delta \rho n1,m1/2$ and $Zn1,m1+1=\u2212Zn1,m1\u22121=\delta \rho n1,m1/2$.

In this work, the VMEC code is run in fixed boundary mode, and a single field period is considered to be the approximation of the cylindrical screw pinch. The major radius of the equilibrium is taken to be $R0,0=100$ m, the unperturbed minor radius to be *a* = 1 m, and the perturbation amplitude to be $\delta \rho 1,2=1\xd710\u22124$ m. This choice of aspect ratio ($R0,0/a=100$) along with a field periodicity of 100 attempts to mimic the cylindrical limit, while the perturbation amplitude places us in a linear limit.

Analytic theory provides a framework against which to compare the VMEC solutions. Linear perturbative ideal MHD theory predicts that a resonant current sheet (the so-called Dirac-*δ* current) should form for the problem considered here. The capturing of such current sheets has been used as a means to verify linear codes. However, for a continuous rotational transform profile, which is presently assumed in linear codes, the existence of flux surfaces near the current sheet is only possible in the limit of vanishing perturbation. Allowing for a jump in the rotational transform, $\Delta \iota $, across the resonant rational ensures that flux surfaces are preserved.^{3} More precisely, the solution of Newcomb's equation, which gives the radial profile of the plasma displacement, $\xi (\rho )$, is consistent with the existence of flux surfaces if the starting equilibrium has a rotational transform profile of the form $\iota (s)=0.6\u22120.2s\xb1\Delta \iota /2$, where the ±sign refers to either side of the resonant surface. As the jump in transform gets smaller, the profile $\xi (\rho )$ becomes steeper and flux surfaces start overlapping when $d\xi /d\rho >1$. In fact, a minimum value for $\Delta \iota $ can be derived and shown to be proportional to the applied resonant perturbation.^{3} Thus, Loizu *et al.* have discovered a new class of solution to Newcomb's equation which preserves the nested flux surface constraint, allowing interpretation of the VMEC results.

In order to gauge the behavior of the code, when not being forced to resolve resonant rational surfaces, a series of equilibria are evaluated with non-resonant boundary perturbations (in fixed boundary mode). In these simulations, the profiles were left unchanged from the cases with resonant boundary harmonics. This simplifies comparison with the resonant cases. In order to make the perturbations non-resonant, the sign of the toroidal perturbed harmonic was changed from *n* = 1 to *n* = −1. This was the simplest and most direct way to achieve a non-resonant perturbation comparable to the resonant one. A series of equilibria with ever increasing radial resolution can then be compared to determine the rate of self-convergence of the code. Here, the error metric is

where *ξ _{ns}* is the radial profile of the perturbed harmonic at a given resolution,

*ξ*is the radial profile of the perturbed harmonic at the highest resolution (radially, 2048 grid points), and the sum is over the radial grid points (

_{highres}*NS*). This gives an appropriate measure of the rate of self-convergence, which can be compared against the expected rate of convergence (for the finite difference numerical method). In the case of VMEC, this is a first order finite difference, suggesting a second order scaling $O(n\u22122)$ (where

*n*is the number of radial grid points).

Self-convergence studies with respect to poloidal and toroidal modes numbers were also conducted. These studies were done using the 512 radial grid point case. To represent our problem, the minimum toroidal mode numbers required are $n=[\u22121,1]$ and poloidal mode numbers $m=[0,3]$. Examination of the resonant response across the rational surface was used as a metric. Maximum toroidal mode numbers from $|n|=[1,4]$ and maximum poloidal mode number $m=[3,14]$ were considered. Little difference was seen between choices of maximal mode number. In the end, $m=[0,9]$ and $n=[\u22123,3]$ were chosen for the spectrum as these were a good tradeoff between a large spectrum and computational speed.

The overall goal of this work is to evaluate the VMEC plasma response at a rational surface using Loizu's solution to the Newcomb equation. Such a comparison may be made with other codes as well, but for now, we focus our attention on the VMEC code.

## III. RESULTS

Simulations of the resonant “screw-pinch-like” configuration indicated a response, in the perturbed harmonics, around the *q* = 2 surface which depended on both the grid resolution and local shear. Figure 2 depicts the profile of the perturbed harmonic across the entire plasma volume, which becomes steep around the *q* = 2 surface. Associated with this steepening is a current density which peaks just inside the *q* = 2 surface. A clear dependence on grid resolution is present in these simulations. Non-resonant perturbations were also examined showing no such response at the rational surface.

As the radial grid resolution increases, the response at the rational surface becomes more localized, as can be seen in Figure 2. The non-axisymmetric current density becomes peaked just inside the q = 2 surface. This was accompanied by a steeping of the gradient in the perturbed quantity (*ξ*) across the q = 2 surface. Self-convergence analysis (Eq. (6)) indicated that the error was scaling as $n\u22121.8$ (Figure 3). This indicates that the code was converging at a rate slightly slower than that can be expected from the radial finite difference ($O(n\u22122)$). Scans of the aspect ratio and perturbation amplitude (at fixed resolution) were also performed (not shown). These scans indicated little dependence of the normalized response (*ξ*) on aspect ratio or perturbation amplitude. Such analysis suggests that so long as the aspect ratio is greater than 6 and boundary perturbation amplitude is less than $1\xd710\u22122$ m, the simulation is in the large aspect ratio and linear limit (respectively).

The non-resonant perturbations that were investigated showed no response at the rational surface, a result which is consistent with linear theory. Self-convergence tests indicated a scaling of $n\u22122.2$, consistent with a predicted scaling of $O(n\u22122)$ for the first order finite difference. It should be noted that the error in the non-resonant case was approximately three orders of magnitude smaller than that of the resonant cases. The largest source of error in the resonant case comes from inside the *q* = 2 surface, while the error is fairly uniform with radius for the non-resonant case. These self-convergence studies seem to suggest that the code is behaving as expected given its numerical treatment.

A scan in the slope of the rotational transform indicated a sensitivity to shear ($d\iota /d\rho $) in the resonant case. The variation in shear was achieved through variation of the slope (*ι*_{1}) of the *ι* profile (assuming a form $\iota (s)=\iota 0+\iota 1s$). Figure 4 depicts the dependence of the resonant response at the *q* = 2 surface as a function of shear. In these cases, as the *ι* profile steepens, the current density at the rational surface localises. This behaviour is consistent with Loizu's solution to Newcomb's equation, where the dependence on the gradient in the response is proportional to the local shear

where *ξ* is the displacement, *ι* is the rotational transform, $\Delta \iota $ is jump in transform associated with the axisymmetric current sheet, and the subscript *r* denotes evaluation at the resonant surface. Loizu's solution predicts that a minimum discontinuity in the rotational transform is required to preserve the flux surfaces across the resonant rational surface. This minimum discontinuity in rotational transform can be written as^{3}

In VMEC, there is no explicit discontinuity in transform; however, using Equation (7), an effective $\Delta \iota $ can be calculated

where the primes indicate derivatives with respect to *ρ* and the subscript *r* implies evaluation at the rational surface. Figure 4 depicts the a-posteriori calculated $\Delta \iota eff$ as a function of the slope in *ξ* at the rational surface. The resolution and shear scan seem to fit the same trend suggesting that this effective $\Delta \iota $ depends on not only the shear but also the radial resolution. The trend also appears to meet the criterion $|\Delta \iota |>|\Delta \iota min|\u223c4\xd710\u22125$ for the existence of flux surfaces (value from Equation (8)). These calculations suggest that the response at the rational surface is not as localized as would be predicted by linear theory.

The linear screw pinch solutions of Loizu were also fit to the VMEC solutions by varying $\Delta \iota $ in the linear model (Figure 5). In this analysis, the $\Delta \iota $ in the linear model was varied until a best fit, of Loizu's solution, to that of VMEC could be found. A perfect fit is not found in any case, but solutions appear to be qualitatively agreeing. In all cases, the effective $\Delta \iota $ found is greater than the $\Delta \iota min$ necessary to ensure the nested flux surface constraint. It should also be noted that this $\Delta \iota $ is larger than that which can be attributed to the finite difference method (larger than radial grid scale differences). Details of the linear model can be found in Loizu *et al*.^{3}

These results indicate that the VMEC code is calculating a plasma response across rational surfaces that is consistent with a perturbative approach that enforces the existence of nested flux surfaces. The plasma response is qualitatively similar for this simple problem despite the lack of an explicit discontinuity in the rotational transform. Calculations of an effective discontinuity provide values which are consistent with the presence of nested flux surfaces across the rational surface. This suggests that while the accuracy of the code may be further improved, the calculated plasma response is consistent with the ideal MHD plasma response. Exact agreement may be obtained if the discontinuous rotational transform profile could be properly implemented in the VMEC code.

## IV. DISUSSION

In this work, a screw-pinch-like equilibrium was examined using the VMEC 3D ideal MHD code. A resonant response was found which was consistent with a linear plasma response assuming the existence of nested flux surfaces everywhere. A test of the non-resonant response indicated the code possessed expected self-convergence properties. Here, the expectation of convergence was based on the radial finite difference used in the code. The consistency with Loizu's solution to Newcomb's equation came despite the lack of an explicit discontinuity in rotational transform. Calculation of an implicit discontinuity in transform suggested values which were consistent with the enforcement of nested flux surfaces in that theory. The implicit transform discontinuity was always much greater than the minimum required to preserve the nested flux surface constraint. Thus, the VMEC plasma response has been verified against the linear ideal MHD equilibrium theory,^{3} in the limit of continuously nested flux surfaces.

The source of the implicit rotational transform discontinuity in VMEC is not clear. In the cases studied, the value of this implicit $\Delta \iota $ was greater than the minimum value necessary to preserve flux surface everywhere. Thus, VMEC is clearly not near the limit for the non-existence of nested flux surface solutions. The implicit $\Delta \iota $ was also much greater than that which could be ascribed to numerical accuracy. It would be assumed that the radial finite difference implies some grid scale $\Delta \iota $. For the shears and grids studied here, a grid scale $\Delta \iota $ would be between $10\u22123$ and $10\u22124$. The studies preformed did show that the response had a dependence on both radial grid resolution and shear. This suggests that there must be a relation between the grid resolution, the shear, and $\Delta \iota $ which limits the $d\xi /d\rho $ achievable in VMEC.

In the work presented here, a relatively simple equilibrium has been investigated. The question still remains to be explored if these results hold for a more experimentally relevant equilibrium. Figure 6 depicts the non-axisymmetric parallel current density for a DIII-D experimental equilibrium computed with VMEC. In this experiment, an *n* = 1 field was applied using the in-vessel coil set (two rows of six coils each). A lobe in the current density appears well correlated with the q = 2 surface, indicating similar phenomena to our screw-pinch problem. Another stronger lobe can be seen around the *q* = 10/3 (20/6) surface, and again at the *q* = 11/3 (22/6) surface. Here, the *n* = 6 resonance has been chosen as the applied field of the I-coils has this fundamental mode. Clearly, the phenomena of shielding currents at rational surfaces are not limited to the screw-pinch problem.

## ACKNOWLEDGMENTS

The authors would like to thank A. Turnbull, N. Ferraro, J.-K. Park, and Z. Wang for their useful discussions on linear code results. The authors would also like to thank the referee for their valuable commentary.

This manuscript is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences and has been authored by Princeton University under Contract No. DE-AC02-09CH11466 with the U.S. Department of Energy. The publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript or allow others to do so, for United States Government purposes.