The extendedMHD NIMROD code [C. R. Sovinec and J. R. King, J. Comput. Phys. 229, 5803 (2010)] is verified against the idealMHD ELITE code [H. R. Wilson et al., Phys. Plasmas 9, 1277 (2002)] on a diverted tokamak discharge. When the NIMROD model complexity is increased incrementally, resistive and firstorder finiteLarmour radius effects are destabilizing and stabilizing, respectively. The full result is compared to local analytic calculations which are found to overpredict both the resistive destabilization and drift stabilization in comparison to the NIMROD computations.
I. INTRODUCTION
Understanding the complex physics of instabilities in the tokamak pedestal region is essential for ensuring that tokamaks achieve their highest performance. Computation plays a critical role in furthering the knowledge of edge dynamics, including advances in the understanding of typeI edgelocalized modes (ELMs).^{1} Nonlinear computation of ELM dynamics remains challenging, but is becoming tractable with modern computing resources.^{2–8} ELM modeling can be roughly broken into three regimes: the early linear regime, the nonlinear regime, and an energy loss regime. Two important aspects of the linear stage are the stability boundary within an operational space, e.g., pedestal current and pressure gradient, and the linear growth rate for each toroidalmode number, i.e., the growthrate or mode spectrum. The dominant instability drive that sets the linearstability boundary is from the peelingballooning mode (PBM)^{9–11} based on the idealMHD model. The linear spectrum is critical to the next phase, the nonlinear regime as it seeds the dynamics. Nonideal effects can greatly impact the growthrate spectrum and are thus important in determining the nonlinear ELM dynamics. In the final regime, highfidelity simulation of the energy loss requires modeling of the fast conductive and convective transport, and accurate coupling to neutrals and the wall. Stateoftheart computations are beginning to push the boundary of current capabilities to modeling that captures accurate dynamics with these additional effects that are relevant during energy loss.^{7,8,12}
As the subsequent nonlinear evolution and energy loss in simulations are greatly impacted by the amplitudes of initial linear perturbations, a large source of uncertainty arises from the linear phase of ELM evolution. Given the difficulty (or impossibility) of measuring the small perturbations that seed an ELM, a reasonable model is to seed it from unstable PBMs that grow from smallwave perturbations. Although the stability boundary is well described by ideal MHD, details of the mode spectrum depend upon specific nonideal effects. In particular, resistive effects may destabilize and firstorder finiteLarmourradius (FLR) effects may stabilize the modes. These latter effects, part of which is drift stabilization, require twofluid modeling with advanced FLR closures. It is important to verify that codes accurately capture these effects within regimes of interest. In this paper, we study these effects on the linear modes with the extendedMHD NIMROD code.^{13,14}
TypeI ELM stability boundaries are dominated by drives described by ideal MHD as the idealMHD operator is stiff, i.e., near the idealMHD stability boundary, small changes in pressure and currentdensity gradients lead to large changes in the growth rate due to the idealMHD terms becoming large. Thus, only small profile changes are required to cross the boundary, and nonideal modifications to the growth rate become secondary when describing the boundary location within operational phase space. This idealMHD stiffness is a general property that impacts the behavior of a wide class of plasmainstability phenomena. For example, in addition to ELM stability, it also describes the dynamics of disruptions caused by coremode activity.^{15,16} However, typeI ELMs are more complicated than coremode disruptions given the broad toroidalmode spectrum and associated nonlinear coupling. Similar to the studies of coremode disruptions, comprehensive understanding of ELM onset requires evolving the equilibrium on the transport time scale through the instability boundary with a verified code that includes nonideal terms to calculate an accurate initial linearmode spectrum.
MHD codes that study ELM instabilities can be placed into two broad categories: (1) extendedMHD codes (e.g., NIMROD and M3DC1^{17}) and (2) reducedMHDbased codes (e.g., BOUT++^{18} and JOREK^{19}). The definition of extended MHD used in this work includes the full idealMHD force operator (see Ref. 20 and references therein). Using the full operator retains the fast compressionalAlfvén waves that enforce force balance and thus are critical to extendedMHD transport calculations with sources and sinks (e.g., Ref. 21). Relative to reduced modeling, the extendedMHD approach also has advantages in avoiding lowβ approximations and a straightforward implementation of the FLR closures as discussed in Section III. However, the inclusion of the full force operator comes at a cost of increased numerical difficulties when solving the equations associated with the fast waves.^{22–26} Verification for all regimes is important, and previous work^{27} verified the idealMHD terms for a nondiverted, peelingballooning unstable test case. As one expects that transport calculations across the ELM stability boundary coupled with an instability model that captures FLR drift stabilization significantly impacts dynamics, emphasis on verifying linear physics of extendedMHD codes thus naturally follows.
In this work, a verification of the full extendedMHD, an initialvalue NIMROD code with the idealMHD ELITE code^{10,11,28} is performed using the Meudas1case comparison.^{29} This is a diverted tokamak case that is peelingballooning unstable. As we later discuss, the inclusion of the xpoint significantly modifies the mode structure, making the comparison in full magnetic geometry an important additional verification. To assess the impact of the twofluid, FLR effects, we compare our results with drift stabilization to a calculation of drift stabilization on the resistive PBM from Ref. 30. Although the approximations made during the analysis of firstorder FLR PBMs in Ref. 30 preclude quantitative comparison, we make a qualitative comparison of our firstorder FLR computations to analytics.
The paper proceeds as follows: In Sec. II, we summarize the benchmark and give an overview of our methodology. The exhaustive details of the methods used in this comparison, in particular, the mechanism whereby a vacuum region is effectively included through the extendedMHD modeling, are previously available in Refs. 27 and 31 and thus our description of this aspect is brief. In Sec. III, the model complexity is gradually increased from the ideallike model through a resistiveMHD model with experimentally relevant profiles to a fully extendedMHD model which includes both parallelclosure and FLRdrift effects. The resulting impact on the growthrate spectrum is analyzed. Comparisons of the computed driftstabilized growth rates are made to analytic treatments of the PBM in Sec. IV. In Sec. V, the density and the temperature are varied at fixed plasma β (thus, the ion gyroradius, ρ_{i}, is modified) in order to ascertain how the FLR effects are impacted. In this study, the bootstrap current is fixed (the zerothorder effect from variations of density modifies the bootstrap current) and the impact of FLR drift effects is isolated. Sec. VI summarizes this work and places it within the broader context of modeling of ELM relaxation events in the tokamak edge with transport effects.
II. VERIFICATION BENCHMARK
We begin with a study of a high resolution, lowersinglenull, JT60Ulike equilibrium (“Meudas1”), which was originally employed in a benchmark of the MARG2D and ELITE codes,^{32} including a close approach to the Xpoint.^{29} This extends the previous benchmarks^{27} of ELITE and NIMROD as it includes diverted magnetic topology and a higher edge safety factor ( $ q 95 = 6.74$, the safety factor at 95% of the normalized poloidal flux) that leads to increased resolution requirements. An idealMHD limit is achieved in NIMROD by using flat density and resistivity profiles inside the last closed flux surface (LCFS) with small resistivity ( $ S = 10 8$), where S is the Lundquist number ( $ S = \tau R / \tau A$), τ_{A} is the Alfvén time ( $ \tau A = R o / v A$), v_{A} is the Alfvén velocity ( $ B / m i n i \mu 0$), τ_{R} is the resistive diffusion time ( $ \tau R = R o 2 \mu 0 / \eta $), $ R o = 2.936 m$ is the radius of the magnetic axis, η is the electrical resistivity, μ_{0} is the permeability of free space, $ m \alpha $ is a species mass (the α subscript denotes ions or electrons in this work), and $ n \alpha $ is a species density. The deuteron mass (m_{i} = 3.34 × 10^{−27} kg) is used. In order to reproduce the vacuum response model outside the LCFS that is used by ELITE, a low density (0.01 of the core density) and a high resistivity (10^{7} times the core resistivity) are prescribed beyond the LCFS (more details on these approximations are in Ref. 27).
The normalized growth rates ( $ \gamma \tau A$ where the linearized mode grows as $ exp [ \gamma t ]$) vs. toroidal mode number ( $ n \varphi $) from NIMROD and ELITE are plotted in Fig. 1. There is good agreement between the codes except for $ n \varphi = 4$ where there is a 27% relative difference. All other cases have a relative difference of less than 8% with typical differences of 5%. The NIMROD convergence in terms of the maximum polynomial order of the spectral elements is shown in Fig. 2. Convergence is most challenging at high wavenumbers where the resolution requirements are most stringent (the poloidal mesh is composed of 72 × 512 spectral elements). These cases converge from the unstable side where the growth rate decreases with enhanced resolution. Thus, the excellent agreement between NIMROD and ELITE at high $ n \varphi $ in Fig. 1 may be spurious and indicate that slightly more resolution is required for $ n \varphi > 25$; however, the shown growth rates are likely within 5% of their converged values. Studying nearly ideal cases with the extended MHD codes such as NIMROD is challenging given the vanishingly small dissipation operators, and convergence is achieved more quickly with the additional nonideal terms in the extendedMHD equations, as in the cases in Sec. III.
Figure 3 shows a poloidal cross section of the magnetic (B_{R}) eigenmode. The mode develops an “interferencepattern” structure near the Xpoint when inboard and outboard fingerlike structures overlap. The finiteelementmesh nodes are superimposed atop the smallestscale subfigure. As established in Fig. 2, this simulation is spatially and temporally converged. The high resolution required to resolve these highq_{95}, diverted cases motivated the development of memoryscaling improvements in the NIMROD code.
III. INCREASING THE MODEL COMPLEXITY
Initialvalue extendedMHD computations, while more computationally expensive than idealMHD calculations with ELITE, are able to include additional effects that are more representative of experimental conditions:

Instead of a magnetostatic vacuum outside the LCFS, NIMROD includes a coldplasma region.

Finitedissipation effects such as resistivity, viscosity, and thermal conduction.

Density and temperature profiles that prescribe the associated drifts and resistivity profile.

Parallelclosure effects that represent large diffusivities oriented along the magnetic field.

Firstorder FLRclosure and twofluid effects.
The authors in Ref. 31 investigated the first three of these effects with the Meudas1 case. We further this work by including the last two effects. The first three effects are also reinvestigated and different modifications to the growth rates are found. These differences are expected as slightly different profiles for density and temperature are used in order to avoid spurious electrostatic modes with the twofluid, firstorder FLR model. The authors in Ref. 31 avoided these modes by employing only the singlefluid resistiveMHD model.
In calculations beyond the ideal model, the density profile is $ n e ( \psi n ) = n e 0 ( p ( \psi n ) / p 0 ) 0.3$ where $ n e 0 = 6 \xd7 10 19 m \u2212 3$ and ψ_{n} is the normalized poloidal flux. Here, p is the pressure as prescribed by the ideal gas law, $ p = n e T e + n i T i$, where $ T \alpha $ is the species' temperature. This choice does not match Ref. 31, but effectively avoids spurious electrostatic modes with the twofluid, firstorder FLR model. These spurious modes are discussed in more detail in Sec. V. The pressureprofile gradient is specified by the Meudas1 case, and setting the edge temperature to 200 eV leads to a core temperature of 5806 eV. The edge temperature is modified independent of the MHDstability (determined by $ p \u2032$) by the transformation $ p ( \psi n ) \u2192 p ( \psi n ) + p c$ where p_{c} is a constant. Without scrapeofflayer (SOL) profiles of density and temperature, which are not included in standard reconstructions, this choice of the edge temperature also specifies the temperature outside the LCFS. The SOL profiles are required to get both the correct profiles outside the LCFS, which determines the vacuum response, and simultaneously set the correct resistivity at the mode resonant location when Spitzer resistivity is used. The nonideal viscous, conductive, and particle diffusivities are set to a small value, one tenth of the resistivity at the LCFS.
Figure 4 plots the growth rates computed by NIMROD on the “Meudas1” case with four different models: ideal MHD (ideal); resistive MHD with reconstructed density, temperature; Spitzerresistivity profiles (1f); and the singlefluid model with parallel Braginskii closures (1fPar) where the parallel viscosity is
and the parallel heatflux vector is
and a twofluid model with parallel closures that includes ion gyroviscosity
where W is the rateofstrain tensor
a fully extended MHD Ohm's law
and separate temperature evolutions with crossheat fluxes
(2fParFLR). Here, v_{i} is the bulkion flow, $ \nu \u2225 i$ is the ion parallel diffusivity, $ \chi \u2225 \alpha $ is the species' parallel diffusivity, $ q \alpha $ (e) is the species' (electron) electric charge, J is the current density, $ b \u0302$ is the magnetic unit vector ( $ b \u0302 = B / B$), and I is the identity tensor. An effective ion charge (Z) of three is assumed. Each model builds on the previously listed and includes all prior terms. Electron viscosity is not included.
When a density profile is included, the normalized growth rate of the mode remains constant with the singlefluid model when the Alfvén time is computed at the mode resonant surface. Our calculations are normalized by the Alfvén time as computed with the values at the magnetic axis. These values remain constant when a density profile is included such that the normalized mode growth rate is effectively increased.
The effect of the resistivity profile is more nuanced. Finite resistivity allows for reconnection within the model by relaxing the frozenflux constraint. This permits resistiveballooning modes that grow faster than their ideal counterparts. Alternatively, the resistive dissipation can stabilize the modes by acting as a dissipative term and modifying the response outside the LCFS from that of a vacuum to that of a plasma. Figure 5 of Ref. 31 shows that the plasma response is stabilizing relative to a vacuum region. Overall, these effects stabilize the ballooning modes at low $ n \varphi $ and destabilize the modes at high $ n \varphi $ as seen in Fig. 4 when comparing the single fluid and ideal normalized growth rates. The effect of including small particle, viscous, and thermal diffusivities on the resistiveMHD mode is small (not shown).
Including the Braginskii parallel closures with large coefficients ( $ \chi \u2225 e = 10 9 \xd7 \eta 0 / \mu 0$ and $ \chi \u2225 i = \nu \u2225 i = 10 8 \xd7 \eta 0 / \mu 0$) has little effect on the growth rates. When thermal conduction is large, the temperature quickly equilibrates along field lines to produce an isothermal response. Thus, NIMROD computations that show little effect from the parallel closures are consistent with the ELITE results that find the shape of the mode growthrate spectrum is largely not sensitive to the value of the ratio of specific heats, Γ, as seen in Fig. 1. When the response to the compressible motion associated with the sound wave is eliminated (Γ = 0 in the ELITE computations), the growth rate is slightly enhanced relative to the adiabatic limit. If the growth rates are modified beyond this small effect, it would be an indication that changes in the energy equation can impact the MHD response, and thus that closure effects are significant. Although the shortmeanfree path Braginskiilike closure^{34,35} is used in our NIMROD computations, we note that other closures, such as longmeanfree path^{36} or the general approach of solving driftkinetic equation^{37} for a closure, can be applied. However, for this case where the parallelclosure effects are negligible, different parallel closures are unlikely to significantly modify the results.
With the full twofluid, FLR model, there is a stabilizing effect on the intermediate and high $ n \varphi $ modes as expected from analytic treatments.^{30} A common reducedMHD approximation for firstorder FLR closures is to assume that the iongyroviscous force and the advection by the diamagnetic drift exactly cancel (this is colloquially known as the gyroviscous cancellation, see, e.g., Ref. 38). However, as discussed by Ramos in Ref. 39, “these cancellations are only partial and not very useful in practice for general magnetic geometries…”. This cancellation is only valid with a large, uniform guide field, without curvature terms (slab approximation), and in the lowβ limit (see, e.g., Ref. 40). The use of the full gyroviscous operator is critical as in addition to the effects that appear as diamagnetic drifts, it also contains terms proportional to firstorder FLR magneticcurvature and gradB drifts.^{41} Our modeling contains not just the firstorder FLR drifts from ion gyroviscosity but also the drifts from the crossheat fluxes that enter the equations on the same order.
Just as the large balance of the idealMHD forces can cause numerical difficulties with the idealMHD force operator, the partial gyroviscous cancellation leads to similar numerical pitfalls and verification is important. The implementation of these terms in NIMROD is verified in cylindrical and slab magnetic geometries against analytic calculations for tearing^{14,40,41} and iontemperaturegradient (ITG)^{42} modes. To test the effect of these terms against analytic theory in full tokamak magnetic geometry, we perform a qualitative comparison with analytic theory in Sec. IV. Full verification is precluded by the approximations made in the analytic theories.
The firstorder FLR model is valid for the toroidal modes presented, although this is not necessarily the case for general edge mode modeling. Assuming $ k \u2243 ( q / r + 1 / R ) n \varphi $ and evaluating local quantities at the peak of the mode eigenfunction ( $ \psi n = 0.969$) on the outboard midplane, the normalized iongyroradius ( $ \rho i = \Gamma m i T i / Z e B$) as a function of toroidal mode is $ k \rho i \u2243 0.0055 n \varphi $. If the firstorder FLR assumption is violated, morecomplex fullionorbit kinetic modeling is required to simulate largefluctuation ELM dynamics. For this case, computations produce approximately the same results with and without the crossheat fluxes. This is consistent with drift analytics^{40} that shows that the crossheat flux is only significant at very low values of plasma β (here the $ \beta = 2 \mu 0 p / B 2$ is 0.1%). Perhaps, coincidentally for this case, the effects of the resistive destabilization and drift stabilization largely counteract one another. The drift stabilized growth rates are not less than the ideal calculation until $ n \varphi \u2273 24$.
IV. COMPARISON TO DRIFT ANALYTICS
To gain insight into our twofluid computations, we compare with analytic descriptions of drift stabilization of the ballooning mode. The dispersion relation for drift stabilizations of the ideal, incompressible ballooning mode is
where γ_{I} is the ideal growth rate in the absence of drift effects and ω is the complex frequency ( $ \omega = i \gamma $).^{43} The diamagnetic drift velocity is typically defined as the perpendiculartothemagneticfield flow contribution from the pressure gradient within the context of a twofluid or generalized Ohm's law. For a tokamak, neoclassical effects damp the poloidal diamagnetic flow contribution. Thus in our case, the remaining toroidal diamagnetic flow determines the diamagnetic frequency, $ \omega * i = k \xb7 v * i = n \varphi / ( n i q i ) \u2202 p i / \u2202 \psi $, where k is the mode wavenumber. Figure 5 shows the comparison to this dispersion relation with full and reduced diamagnetic drift values. The value of $ \omega * i$ is computed at the radial location where the mode structure peaks. The normalized frequency from the diamagnetic drift velocity, $ \omega * i \tau A$, varies linearly with toroidal mode number as $ 0.00483 n \varphi $. This “local” approximation, which neglects the radial variation of the $ \omega * i$ profile, is known to overestimate the effect of drift stabilization.^{44,45} We find that a reduction of the value of $ \omega * i$ (approximately by a factor of 2–4) gives the best agreement between Eq. (7) and the twofluid NIMROD calculation. A similar comparison of ELITE to the twofluid BOUT++ model on a different case finds the drift stabilization is overpredicted^{45} by a comparable factor.
One of the limitations of Eq. (7) is that it does not include the effect of resistive destabilization on the high $ n \varphi $ modes, as is included in our twofluid computations. The analytic treatment of Ref. 30 provides a dispersion relation (Eq. (43) of the reference; referred to as HRP Eq. (43) in this discussion) that includes the effect of resistive destabilization in addition to drift effects. In order to compare with this theoretical treatment, parameters are computed on the outboard midplane at the radial location where the mode structure peaks, consistent with the previous $ \omega * i$ calculation. This yields local values of the safety factor, q = 7.34, normalized magnetic shear, s = 8.8, and normalized pressure gradient, $ \alpha = 28.8$, where the definitions of Ref. 30 are used. The $ \Delta B \u2032$ stability parameter is determined by using the ideal growth rate ( $ \Gamma = 5 / 3$) and solving Eq. 35 of Ref. 30. $ \Delta B \u2032$ varies with toroidal mode number but falls in the range of −17.4 to −12.1.
Comparison of HRP Eq. (43) with $ \omega * = 0$ and the NIMROD resistiveMHD results differ by a factor approximately 2.5 at large toroidal mode numbers where the analytic calculation results in greater growth rates and there is no stabilization at low $ n \varphi $. As the analytic calculation is a local calculation, it does not capture the stabilizing influence on the global mode of finite resistivity, in particular, it misses the effect of modifying the response outside the LCFS from that of a vacuum to that of a resistive plasma. Our interest is largely to compare and contrast the predicted drift stabilization, and thus, we reduce the resistivity to 35% of the value used in the NIMROD computations when calculating the analytic values of HRP Eq. (43). This produces growth rates with $ \omega * = 0$ that roughly match the NIMROD resistiveMHD computations at high $ n \varphi $ as shown in Fig. 6.
More significantly, Fig. 6 compares the results of our linear NIMROD computations with the driftstabilized growth rates computed from Eq. (43) of Ref. 30. Again, the analytics of HRP Eq. (43) overestimate the effect of drift stabilization where the reduction of $ \omega *$ that leads HRP, Eq. (43), to qualitatively agree with the twofluid NIMROD computations is between 2 and 4 (dependent on $ n \varphi $). There is little effect on the low $ n \varphi $ modes and the ion driftwave resonance predicted in Ref. 30 is not observed.
V. FLR PARAMETER SCAN
In order to ascertain the role of twofluid effects further, we vary the density and temperature while scaling parameters to maintain constant β and S (other dissipation parameters are scaled relative to resistivity). In experimental discharges, the zerothorder effect on the edge of modifying the density is to modify the plasma collisionality and thus the bootstrap current. This causes a transition from ballooninglike modes at high density to peelinglike modes at low density. Our computations use a fixed equilibrium current and thus are not sensitive to this effect. Instead, we isolate the effect of the modification of ρ_{i} (and the associated ion skin depth, d_{i}) within the context of our twofluid, firstorder FLR model. This contrasts with Ref. 46 that examines these effects in concert and Ref. 47 that considers only the effect of modifications to the current profile.
Figure 7 shows the linear growth rate from NIMROD computations for resistiveMHD and three twofluid cases with varying densities (the core density is varied from $ 3 \xd7 10 19 m \u2212 3$ to $ 1.2 \xd7 10 20 m \u2212 3$ where $ 6 \xd7 10 19 m \u2212 3$ is the value used in the prior computations discussed in the paper). The effect of driftstabilization is stronger at low densities (ρ_{i} scales as $ 1 / n i$ at constant β) as expected. The cases with $ n e = 3 \xd7 10 19 m \u2212 3$ are not plotted for $ n \varphi > 18$ as the most unstable mode is no longer the peelingballooning mode but rather a dominantly electrostatic mode related to the ITG as studied in the context NIMROD twofluid advance in Ref. 42. This region is excluded as examination of ITG mode dynamics is outside the scope of this work.
VI. DISCUSSION AND CONCLUSIONS
Our comparisons to the ideal PBM growth rates computed by the ELITE code and analytic descriptions of the resistive, twofluid PBM provides another verification test of the NIMROD code as a means to study the tokamak edge. Quantitative agreement with ELITE is achieved for these cases, extending prior work^{27} that did not include diverted magnetic topology.
Twofluid and FLRdrift effects can both produce stabilizing (drift stabilization through the sheared motion of the ion and electron fluid responses) and destabilizing (species decoupling and enhanced growth through mediation by the more mobile electron fluid) effects (see, e.g., Refs. 40 and 41). At intermediate and large mode numbers, the inclusion of finiteresistivity is destabilizing, while the FLRdrift effects are stabilizing. These effects are largely offsetting, resulting in a toroidalmode growthrate spectrum with the full extendedMHD model that resembles the ideal calculation. This result is not general but rather case dependent. Consistent with ELITE calculations that vary the ratio of specific heats from the adiabatic through the isothermal limit, we find the addition of anisotropic thermal conduction produces only a small modification of the growth rate.
Relative to the NIMROD calculations, analytics^{30} overpredict both the effects of resistive destabilization and FLRdrift stabilization for the case studied. Qualitative agreement with the analytic descriptions of the peelingballooning mode is found where high $ n \varphi $ modes are susceptible to resistive destabilization and drift stabilization in both the analytics and NIMROD computations. One of the limits of the analytic derivation is a “local” approximation which leads to the over prediction of FLRdrift stabilization.^{44} Given the highly complex nature of both the configuration of the tokamak edge and the twofluid model equations, quantitative agreement is not expected as the local approximation neglects the global configuration and profile effects.
Regarding our investigation on the impact of the size of the ion gyroradius within the context of NIMROD's full firstorder FLR, twofluid model, we find significant changes on the growth rates of high $ n \varphi $ modes at varying densities but constant β. However, the most unstable mode number is largely unchanged. Thus, we conclude the zerothorder effect of modifications to the bootstrap current largely dominates the determination of the most unstable mode during collisionality scans, while lower densities (large ρ_{i}) increase the magnitude of the drift stabilization on the high $ n \varphi $ modes consistent with the discussion of Ref. 48.
ACKNOWLEDGMENTS
We thank Carl Sovinec, Chris Hegna, and Nate Ferraro for discussions involving this paper and N. Aiba for providing the Meudas1 equilibrium. This material is based on the work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Award Nos. DEFC0206ER54875 and DEFG0208ER54972. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DEAC0205CH11231.