A Parametric Study of Extended-MHD Drift Tearing
aa r X i v : . [ phy s i c s . p l a s m - ph ] J u l A Parametric Study of Extended-MHD Drift Tearing
J. R. King and S. E. Kruger Tech-X Corporation, 5621 Arapahoe Ave. Suite A Boulder, CO 80303 (Dated: draft August 22, 2018)The linear drift-tearing mode is analyzed for different regimes of the plasma- β , ion-skin-depthparameter space with an unreduced, extended-MHD model. New dispersion relations are found atmoderate plasma β and previous drift-tearing results are classified as applicable at small plasma β .The drift stabilization of the mode in the regimes varies from non-existent/weak to complete. As thediamagnetic-drift frequency is proportional to the plasma β , verification exercises with unreduced,extended-MHD models in the small plasma- β regimes are impractical. The new dispersion relationsin the moderate plasma- β regimes are used to verify the extended-MHD implementation of theNIMROD code [C. R. Sovinec et al ., J. Comput. Phys. , 355 (2004)]. Given the small boundary-layer skin depth, discussion of the validity of the first-order finite-Larmour-radius model is presented. PACS numbers: 52.30.-q, 52.65.-y, 52.35.-g, 52.55.Tn, 52.40.Hf, 02.60.Lj, 52.35.VdKeywords: Drift-tearing mode, Extended-MHD modeling, Code Verification
I. INTRODUCTION
Experimental, fusion-plasma discharges typically operate in regimes away from ideal-magnetohydrodynamic(MHD) stability boundaries. The ideal-MHD modes that exist outside these boundaries, which are unableto modify the magnetic topology, are often deleterious to confinement and can lead to a rapid loss of theplasma stored energy. Analysis with a resistive-MHD model shows a second class of modes are possible.These resistive-MHD modes are a combination of macroscopic ideal-MHD behavior through-out most of theplasma volume and boundary-layer dynamics where resistivity is important near a resonant magnetic-fluxsurface, a surface where the mode structure and the magnetic topology are aligned in poloidal and toroidalperiodic variation. Although the plasma dynamics associated with these resistive modes are usually lessviolent than ideal modes, finite resistivity allows for modification of the magnetic topology. For example,magnetic islands formed from saturated resistive-tearing modes can enhance energy and particle transportfrom the plasma core to the edge via large field-aligned transport.The tearing instability [1] is one such multi-scale mode: a combination of macroscopic structure, theideal-MHD response through-out most of the plasma volume; and microscopic structure, the boundary-layerphysics near the resonant surface which minimally includes resistive MHD. Ideal-MHD flows advect magneticflux to the resonant surface where a large, localized current sheet is formed. This leads to slow growth ona hybrid-time scale that is a combination of the ideal Alfvén time and the time scale of the pertinentboundary-layer physics. With a resistive-MHD model, the current-sheet size is determined by the magnitudeof the plasma resistivity: smaller resistivity results in a more localized layer. In high-temperature fusionplasmas, which have very small resistivity, the boundary-layer width can approach the ion gyroradius wherefinite-Larmour-radius (FLR) and electron-ion-fluid-decoupling effects become important. When the moremobile electron fluid is decoupled from the ion fluid near the layer, it can more effectively transport fluxinto the layer and thus destabilize the mode (increase the growth rate). Alternatively, when the fluids aredecoupled and drift in opposed directions within the resonant flux surface, the sheared relative motion canstabilize the mode (reduce the growth rate). A sufficient model to capture these FLR effects to first orderis extended-MHD with Braginskii-like closures [2–5]. The zeroth-order plasma drift, the E × B drift, causesthe electron and ion fluids to drift with the same velocity and thus are not stabilizing. The first-order FLRdrifts have a orientation that is dependent on the sign of the charge of the species and thus are stabilizing.With respect to influence on the tearing mode, the most studied first-order FLR drift is the fluid diamagneticdrift [6] but stabilizing effects are also attributed to drifts proportional to the gradient and curvature of themagnetic field [7].A previous parametric regime analysis of the tearing mode without drift effects is given by Ahedo andRamos [8]. They characterize small- ∆ ′ tearing-mode parameter space by seven regimes as illustrated Table 1 parameter spacePR1 (single (cid:1) uid)PR0 (no solution)PR2 PR6PR3 (electron MHD)PR4
Figure 1: Tearing mode parameter space in terms of normalized β , ¯ τ − Q , and d i , ¯ σ , (as originally defined in Ref. [8]).Growth rates from PR1 through PR5 are used within the normalizations as appropriate. The colored (dark) areamaps the region of interest for tokamak fusion plasmas with parameter ranges as defined in Tab. I. A first-orderion-FLR model is valid in the blue, dotted region ( ρ i < . δ ), and invalid in the solid, red region ( ρ i > . δ ). Thereis a wavy, purple region which contains both valid and invalid cases as the normalized parameter space and ρ i /δ donot have a one-to-one mapping. The diagram uses a small parameter value of . to determine regime boundaries.Points A through D correspond to the ω ∗ → limit of the verification scans of Sec. VI and the modification of ¯ τ and ¯ σ as ω ∗ is increased is illustrated with dashed lines for the ranges of ω ∗ included in the verification exercises. schematically in Fig. 1. A single-fluid model (resistive MHD) describes the dynamics in parameter-space-region PR1 as first discovered by Furth et al. [1]. In PR5, at very small values of β and large ion-skin depth( d i ), the semicollisional description of Drake and Lee is valid [9]. Without drift-effects, the semicollisionaldescription is valid when β is smaller than the square of the tearing skin depth normalized by the modewavenumber (ignoring some factors of order unity). Thus even for a large tearing skin depth of 1cm, thevalidity constraint is still approximately β << − with a mode wavelength of 1m and it is unlikely thisregime is relevant to tokamak discharges. At moderate values of the plasma- β parameter and d i (thus mod-erate values of the ion gyroradius, ρ i ∼ √ βd i ) the tearing dispersion relation from electron-MHD [10] isrecovered (PR3). Mirnov et al. derive a unified dispersion relation for PR4 which limits to that found inPR3 and PR5 [11]. They describe the decoupling effects of the mode as mediated through interaction withthe kinetic-Alfvén wave in PR3 and the whistler wave in PR5. The dispersion relations for the remainingtransitional regimes in this parameter space (PR2 and PR6) are derived by Ahedo and Ramos [8]. There isno known solution in PR0.Our results largely follow the parameter-space characterization of Ref. [8], however our calculations includediamagnetic and magnetic-field-gradient drift contributions. In a sense, the main concept of our study is toadd a third dimension out of the page of Fig. 1 that corresponds to the drift frequency. In Sec. II, we describethe extended-MHD model and our small- ∆ ′ , large-guide-field assumptions. These equations are linearizedand reduced to a system of two second-order equations in Secs. III and IV. Our intention is to clarify therelevant regimes to fusion plasmas and benchmark extended-MHD drift-tearing computations which use anunreduced-MHD model. As such, our study differs from much of the prior work in that we do not startwith a reduced-MHD model, but rather we apply tearing ordering to the full extended-MHD equations. Ourmain dispersion relation results are derived in Sec. V for drift tearing in PR1 through PR5. We recover theresult of Coppi at small values of the plasma- β parameter in the single-fluid regime (PR1, Refs. [6, 12]) andthe result of Drake and Lee in the semicollisional regime (PR5, Ref. [9]). New dispersion relations are foundin PR2 through PR4.The linear drift-tearing response in the moderate- β regimes is necessary to verify extended-MHD codes atthe parameters of typical use cases. In Sec. VI we present the results of a verification exercise between theNIMROD extended-MHD code [13] and our new drift-tearing dispersion relations in PR2 through PR4. The parameter range S - k ⊥ d i . - . β . - . k − ⊥ ∆ ′ ǫ B S , the ion skin depth, d i , the plasma β , the tearing stability parameter, ∆ ′ , and the guide-to-sheared magnetic-field ratio, ǫ B (definitions are provided in the text). extended-MHD drift terms are complicated even in primitive form, as described in Sec. II. The complexityof the model makes implementation and numerical analysis challenging and thus makes verification all themore important. The NIMROD extended-MHD algorithm has previously been benchmarked in PR1 throughPR5 without drift effects [14]. Our verification exercise with drift effects extends this benchmark to moreexperimentally relevant regimes.Our new verification scans of the drift frequency begin at points A-D in the ¯ τ Q - ¯ σ parameter space of Fig. 1.As the drift frequency increases, the ¯ τ Q parameter of Fig. 1 increases and the single-fluid tearing skin depth, δ , decreases. Thus cases move within the ¯ τ Q - ¯ σ parameter space of the figure down and slightly to the rightas illustrated in the figure by the dashed lines. As the diamagnetic-drift frequency scales proportionally to β , we have found it difficult to satisfy all of the analytic asymptotic-limit requirements of the tearing modewhile running with appreciable drift frequency for a verification exercise purely in the low- β semicollisionalregime (PR5).If the effect of electron inertia is larger than that of resistivity, the dynamics are described by so-calledcollisionless physics. In this regime, the growth rate is independent of resistivity. Electron inertia scalesproportionally to the electron skin depth squared, where the ion and electron skin depths ( d i and d e ) havea fixed ratio equal to the square root of the mass ratio. Thus in PR3 through PR5 as d i is increased inFig. 1, the mode will ultimately become collisionless. Also for this reason it is difficult, if not impossible,to compose collisionless cases in PR6, PR1 and PR2 unless one is using a model with an enhanced electronmass or operating with extremely low plasma β . We interpret our drift-tearing results in the transition-to-collisionless electron-MHD regime (PR3) and discuss implications for extended-MHD modeling in Sec. VI.Fitzpatrick points out that when one compares the size of the electron gyroradius, ρ e , to the single-fluidtearing skin depth, δ , in the moderate- β regimes (PR3-5) for collisionless cases, a first-order electron-FLRmodel is invalid throughout all of the electron-MHD regime (PR3) and much of PR4 as ρ e > δ . However themodel is valid in the semicollisional regime (PR5) [15]. A corollary to this argument is that a first-order ion-FLR model will be invalid when ρ i > δ . Consider the expected fusion-plasma parameters as listed in TableI; the limits in the ¯ τ Q - ¯ σ parameter space defined by the Tab. I parameters are superimposed onto Fig. 1. Afirst-order ion-FLR model is valid when √ βd i ∼ ρ i << δ , which encompasses many of the fusion-relevantcases in PR0, PR1 and PR2. For these parameters, a first-order electron-FLR model is always valid (witha realistic mass ratio, µ = m e /m i ) as ρ e /δ = √ µρ i /δ ∼ √ µβ ¯ σ and ¯ σ never exceeds a value of 100 with theparameters of Tab. I. We note two reasons for studying drift tearing with a first-order ion-FLR model outsideits regime of strict validity. First, the model may be outside the region of strict validity only for linear modes.With nonlinear dynamics, the tearing skin depth is no longer a well defined concept and, strictly from thelinear definition, it broadens as the mode approaches saturation. In these nonlinear regimes, model validityis determined largely by the constraint kρ i << that is more easily satisfied by long-wavelength tearinginstabilities ( k is the perturbation wavenumber). First-order FLR, extended-MHD modeling is typicallyinterested in the nonlinear evolution of the plasma; however, most computations first encounter a lineargrowth phase that is still important to both understand and ensure that it is calculated correctly. Second,the mode dynamics transition to an electron-MHD description as ρ i becomes large and the ion fluid becomesdemagnetized on the small tearing-skin-depth scale and decoupled from the electron fluid. If a first-orderFLR model is capable of correctly modeling these electron-fluid dynamics, it may be qualitatively descriptiveof the electron dynamics outside its regime of strict validity. Qualitatively descriptive but computationallytractable first-order FLR, extended-MHD modeling is preferable to modeling with full-orbit ion dynamicswhen the latter is computational intractable. II. MODEL EQUATIONS AND ORDERINGS
With an unreduced-MHD model, the plasma fluid is described by a continuity equation, ∂n∂t = −∇ · n v , (1)for the plasma density ( n ) evolution, a center-of-mass momentum equation, m i n d v dt = J × B − ∇ p − ∇ · Π i , (2)for the bulk-plasma velocity ( v ), and an energy equation, n Γ − d α T α dt = − p α ∇ · v α − ∇ · q α , (3)for the plasma temperature ( T α ). The subscript indicates either the ion or electron species, m α is a species’mass, and Γ is the adiabatic index. The plasma is assumed to be an ideal gas and thus the species pressure( p α ; p = P p α ) is given by the ideal-gas law, p α = nT α . As appropriate for low-frequency plasma dynamics,we assume quasi-neutrality ( n e ≃ n i for an ion charge state of unity) and drop the displacement-current termin Ampere’s law ( µ J = ∇ × B where µ is permeability of free space), which provides a relation betweenthe magnetic field ( B ) and the current density ( J = ne ( v i − v e ) where e is the electron charge). Theseapproximations analytically eliminate both light and Langmuir waves. The electron momentum equation isused as an expression for the electric field ( E ), E = − v × B + J × B ne − ∇ p e ne − ∇ · Π e ne + η J − m e e d e v e dt , (4)commonly referred to as the generalized Ohm’s law ( m e is the electron mass and η is the electrical resistivitycaused by electron-ion collisions). Faraday’s law ( ∂ B /∂t = −∇ × E ) in conjunction with Eqn. (4) producesthe induction equation, which describes the evolution of the magnetic field. This system of equations isconsidered to be a two-fluid model when the Hall term ( J × B /ne ) is retained as the magnetic field is thenadvected by the electron flow ( v e = v i − J /ne ) instead of bulk-flow advection from the v × B term.These equations require closure expressions for the stress tensors ( Π α ) and heat fluxes ( q α ). We use theBraginskii-like [3–5] ‘cross’ terms (first-order FLR terms) as the closure: gyroviscosity, Π α = m α p α q α B h ˆ b × W α · (cid:16) I + 3ˆ b ˆ b (cid:17) − (cid:16) I + 3ˆ b ˆ b (cid:17) · W α × ˆ b i , (5)and cross-heat flux, q = 5 p α q α B ˆ b × ∇ T α , (6)where q α is a species’ charge. The rate-of-strain tensor ( W α ) is defined as W α = ∇ v α + ∇ v Tα − (2 / I ∇ · v α . This choice of closure neglects the perpendicular and parallel (to B ) closure terms and additionalcontributions to the gyroviscous stress [3, 16]; however, the retained terms are commonly included in state-of-the-art extended-MHD codes and have contributions that enter the model equations on the same order asthe diamagnetic-drift terms.To further estimate the importance of the cross-closure terms, consider flows on the order of the soundspeed, c s = p Γ ( T i + T e ) /m i , which for comparable species’ temperatures is on the same order as theion thermal speed, v T α = p T α /m α . The ion gyroviscous term then scales as ρ i /L relative to the ∇ p term in the momentum equation, Eqn. (2), whereas the electron gyroviscous term scales as p m e /m i ( ρ e /L ) relative to the ∇ p e term in the generalized Ohm’s law, Eqn. (4). Here ρ α = v T α /ω cα is the gyroradius where ω cα = q α B/m α is the gyrofrequency and L is a characteristic gradient length scale. Furthermore, the ratio ofthe electron to ion gyroradius is the square root of the mass ratio, p m e /m i . Thus if the ion gyroviscous termis significant and the first-order ion-FLR model remains valid, ρ i /L . O (1) , then the electron gyroviscousterm is expected to be smaller than other terms in the generalized Ohm’s law by at least the mass ratio.As such, we neglect contributions from electron gyroviscosity in our equations. This assumption leads to abreak-down of the model for the collisionless drift-tearing mode where these scalings do not apply withinthe layer, as discussed in Sec. VI. Next consider the cross-heat flux terms relative to p α ∇ · v α in Eqn. (3).The ion cross-heat flux scales as ρ i /L , but the electron cross-heat flux scales as p m i /m e ( ρ e /L ) . Thus ifthe ion cross-heat flux is significant ( ρ i /L . O (1) ) then the electron cross-heat flux enters the equations onthe same order and must be retained.For the purposes of our study, the tearing instability is generated from an imposed ˆ z -oriented current sheetin a Cartesian slab. There are distant conducting walls at x = ±∞ , the ˆ y and ˆ z directions are infinite, andthe ˆ z direction is symmetric. The tearing mode drive is fueled by free energy from the global configurationbut growth the of the mode is limited by the small-scale physics that breaks the frozen-flux theorem withinthe tearing boundary layer. As this boundary-layer physics is the focus of our study, the slab configurationis locally analogous to a toroidal configuration without curvature contributions where ˆ x is a radial (flux)coordinate, ˆ y is approximately a cross-field coordinate and ˆ z is approximately a parallel-field coordinate. Wedecompose all fields into imposed, x -dependent, background fields (‘0’ subscript) and periodic-in- ˆ y , pertur-bation fields (tilde), e.g. B = B ( x ) + ˜ B ( x ) exp ( iky + γt ) . Here k = k ˆ y is the perturbation wavenumberand γ is the complex growth rate. The radial ( ˆ x ) component of all background vector fields is zero. Pertur-bation vector fields and wavenumber use a magnetic-coordinate system where the ˆ y and ˆ z components areexpressed as parallel-to and perpendicular-to the magnetic field.In our subsequent analysis, we ignore the effects of flow shear but retain the effect of advection by bulkbackground flows. We impose orderings appropriate for the tearing boundary layer: (1) the equilibriummagnetic-shear-length scale ( L s ) is comparable to the inverse wavelength, kL s ∼ O (1) ; (2) a moderatelylarge guide- to shear-magnetic-field ratio, such that ǫ B = B z ( x = 0) /B y ( x = ∞ ) ∼ O (cid:0) ǫ / (cid:1) ; (3) a smalltearing skin depth ( δ ), kx ∼ kδ ∼ O ( ǫ ) ; (4) slow dynamics, ωτ A ∼ O (cid:0) ǫ / (cid:1) ; and (5) slowly varying profileswithin the layer, e.g. δn ′ /n ∼ O ( ǫ ) . Here ǫ is a small parameter ( ǫ << ) and τ − A = kv A = kB / √ m i n µ is the Alfvén time.These assumptions are consistent with the expected conditions for core tearing in a high-temperaturetokamak discharge. Our analysis can accommodate very small values of β ( c s /v A ), however we assume thegrowth rate is subsonic ( γ << k c s ). Ahedo and Ramos show that when this assumption is violated withoutdrift effects, the eigenfunction structure is modified but the growth rate is unchanged [17]. We assume thatthe electron-inertia term is dominated by the contribution from current density and that advection in theelectron inertia term, which is of the same order as electron gyroviscosity, is small. Thus after linearization, m e e d e v e dt ≃ m e e (cid:18) v e · ∇ + ∂∂t (cid:19) ˜ v e ≃ m e n e e ∂ ˜ J ∂t . (7)As the linearized contributions from both electron inertia and resistivity are proportional to ˜ J , we simplifythe subsequent equations by combining these terms and forming a generalized resistivity, η g µ = ηµ + d e γ . (8)The relative magnitude of resistivity compared to electron inertia classifies the tearing mode as collisionless( d e γ >> η/µ ) or collisional ( d e γ << η/µ ) where d α is a species’ skin depth ( p m α /µ ne ). Similarly,we use a generalized Lundquist number, S g = v A µ /k ⊥ η g . In the following discussion, we use two normal-izations: the hat which indicates normalization by Alfvén time/velocity and characteristic field strengths( ˆ ω = ωτ A , ˆ L = kL , ˆ v = v/v A ( x = 0) , ˆ B = B/B ( x = 0) , ˆ n = n/n ( x = 0) , and ˆ p = p/v A m i n ( x = 0) ) andthe overbar which is a tearing specific normalization introduced in Sec. IV. III. LINEARIZED EQUATIONS
Following convention, we define ˆ ξ = ˆ γ ˆ v x as the displacement vector and ˆ Q = ˆ k ˆ B k − i ˆ k k ˆ B ′ x + i ˆ λ ˆ B x , (9)consistent with Ref. [8] where λ = µ J · B /B and ˆ k k = k · B /kB . After linearization and applying theassumptions of Sec. II, the radial induction equation becomes ˆ γ e ˆ B x = i ˆ k k ˆ γ ˆ ξ + ˆ k k ˆ d i ˆ Q + S − g ˆ B ′′ x . (10)The left side of this equation is a term representing the rate-of-change of ˆ B x . The notation ˆ γ i = ˆ γ + i ˆ k · ˆv ,and ˆ γ e = ˆ γ + i ˆ k · ˆv e ≃ ˆ γ i − i ˆ ω ∗ gathers the advective and temporal-derivative contributions into a singleterm. The terms on the right side of Eqn. (10) result from the v × B , Hall, and resistive/inertial terms,respectively. Contributions from the ∇ p e term vanish. Other than ignoring flow shear and applying ourordering to resistive/inertial term, Eqn. (10) is exact.The location where k · B = 0 is the resonant magnetic-flux surface. Away from the resonant surface thecontribution from the v × B term dominates and all other terms may be neglected. When fluid decouplingand/or drift effects are significant, the Hall term dominates near the resonant surface. At the resonantsurface the v × B and Hall terms vanish and thus the resistive and inertial contributions must be retained.Our calculations assume the resonant surface is located at x = 0 . The standard treatment of these equationsis to apply a boundary-layer analysis, where the ideal-MHD equations describe the solution in the outerregion (away from the resonant surface), and the full model is used in the inner layer near the resonantsurface. These solutions are matched using the discontinuity in the logarithmic derivative of the perturbedradial magnetic field of the outer solution ( ∆ ′ ), ∆ ′ = lim ǫ → ˜ B ′ x ( x ) (cid:12)(cid:12)(cid:12) x = ǫx = − ǫ ˜ B x (0) , (11)where the prime indicates a partial derivative with respect to x . With a resistive-MHD model, an equilibriumis tearing unstable ( γ > ) if ∆ ′ > [1]; thus ∆ ′ is both a matching and stability parameter. We assumethat ∆ ′ δ ∼ O (1) and thus ˆ B ′ x ∼ ˆ B x , as follows from Eqn. (11). Expanding ˆ B x at x = 0 , ˆ B x = ˆ B x (0) + ˆ B ′ x (0) ˆ x + ... , (12)and noting that ˆ x ∼ O ( ǫ ) allows us to treat ˆ B x as a constant - an assumption known as the constant- ψ approximation. Derivatives of other perturbed fields are assumed to raise the relative size of the field by ǫ − ,e.g. ǫ ˆ ξ ′′ ∼ ˆ ξ and ǫ ˆ B ′′ x ∼ ˆ B x . This approximation results from the large, localized gradients of perturbedfields within the boundary layer. Consider, for example, that the reconnecting inflows of the tearing modeproduce a displacement vector that changes sign across the boundary layer.After linearization, the parallel induction equation becomes ˆ γ e ˆ B k = − ˆ ∇ ⊥ · ˆ v + ˆ ω ∗ ˆ γ ˆ ξ ˆ d i − (cid:18) i ˆ ω ∗ + i ˆ ω ∗ n Γˆ c s (cid:19) ˆ Q + ˆ k k ˆ d i h ˆ B ′′ x + i ˆ k k ˆ B ′k − ˆ B x i − i ˆ ω ∗ i ˆ n − i ˆ ω ∗ n Γˆ c s ˆ p e + S − g ˆ B ′′k . (13)where ω ∗ α is a species diamagnetic-drift frequency ( kp ′ α /n eB ), ω ∗ is the total diamagnetic-drift frequency( ω ∗ i + ω ∗ e ), ω ∗ n is the density-gradient drift ( kT n ′ /n eB ) and ∇ ⊥ = ∇ − ik k ˆ b . The first two pairs of termson the right side are the contributions from the v × B and Hall terms, respectively. The terms involving ˆ n and ˆ p e result from the ∇ p e term and the last term is the effect of resistivity and electron inertia.The components of the linearized momentum equation are ˆ γ i ˆ γ ˆ ξ = ˆ ω ∗ ˆ d i ˆ B k + i ˆ k k ˆ B x − ˆ B ′k − ˆ p ′ − (cid:16) ˆ ∇ · ˆ Π (cid:17) x , (14) ˆ γ i ˆ v ⊥ = − i ˆ Q − i ˆ p − (cid:16) ˆ ∇ · ˆ Π (cid:17) ⊥ , (15)and ˆ γ i ˆ v k = − ˆ ω ∗ ˆ d i ˆ B x − i ˆ k k ˆ p − (cid:16) ˆ ∇ · ˆ Π (cid:17) k . (16)The perpendicular and parallel components (Eqns. (15) and (16)) are used to construct an expression for ˆ ∇ · ˆ v . The first terms on the right side of Eqns. (14) and (16) are drift contributions from J × B .The linearized continuity, ion-energy and electron-energy equations are ˆ γ i ˆ n = − ˆ ω ∗ n Γˆ c s ˆ γ ˆ ξ ˆ d i − ˆ ∇ · ˆ v , (17) ˆ γ i ˆ p i = − ˆ ω ∗ i ˆ γ ˆ ξ ˆ d i − ˆ c si ˆ ∇ · ˆ v − (Γ −
1) ˆ ∇ · ˆ q i , (18)and ˆ γ pe ˆ p e = − ˆ ω ∗ e ˆ γ ˆ ξ ˆ d i − ˆ c se ˆ ∇ · ˆ v + σ pe ( iω ∗ e − Γ f T e iω ∗ n ) (cid:16) ˆ Q − i ˆ λ ˆ B x (cid:17) − σ pe ˆ c se (cid:16) i ˆ ω ∗ + i ˆ k k ˆ λ ˆ d i (cid:17) ˆ n − (Γ −
1) ˆ ∇ · ˆ q e , (19)respectively. Advection by fast, parallel, electron flows can be computationally expensive to model inextended-MHD computations. A common computational practice is to use the bulk flow in the advec-tive term of the electron-energy equation which circumvents the large computational cost of the fast electronflows. To allow for a systematic study of the effect of different advective models, we introduce the σ pe and ˆ γ pe notation. If the advective term uses the bulk flow then ˆ γ pe = ˆ γ i and σ pe = 0 , whereas advection by theelectron flow leads to ˆ γ pe = ˆ γ e and σ pe = 1 . To compute the linearized cross heat-flux contributions we firstexpand the heat-flux vector as ∇ · q α = ∇ · (cid:20) p α q α B ˆ b × ∇ T α (cid:21) = 5 p α nq α B (cid:20) µ J · (cid:18) p α ∇ nn − ∇ p α (cid:19) − λ B · (cid:18) p α ∇ nn − ∇ p α (cid:19)(cid:21) + 5 p α nq α B h p α n ( B × ∇ n ) − B × ∇ p α i · B · ∇ B + 5 p α n q α B B · ( ∇ p α × ∇ n ) . (20)Noting that J ·∇ f , B ·∇ f , B ·∇ B , and ∇ f ×∇ g vanish for our slab configuration, we may assume thecoefficients of these terms are equilibrium quantities during linearization. After linearization and ordering(specifically, we drop terms where ˆ ω ∗ >> ˆ k k ˆ λ ˆ d i ), we find (Γ −
1) ˆ ∇ · ˆ q α = i ˆ ω ∗ q α ˆ c sα Γ ˆ n − i ˆ ω ∗ qα ˆ p α − (ˆ γ α − i ˆ ω ∗ qα ) ˆ c sα C qα (cid:16) ˆ Q − i ˆ λ ˆ B x + 2 i ˆ k k ˆ B ′ x (cid:17) , (21)where C qα = σ qα i ˆ ω ∗ α − f T α i ˆ ω ∗ n ˆ γ α − i ˆ ω ∗ qα , (22) i ˆ ω ∗ q α = σ qα (cid:0) Γ i ˆ ω ∗ α + ˆ c sα i ˆ ω ∗ (cid:1) , (23)and i ˆ ω ∗ qα = σ qα f T α (cid:0) Γ i ˆ ω ∗ n + ˆ c s i ˆ ω ∗ (cid:1) . (24)Again we introduce σ qα as a marker with value σ qi = − σ qe = (5 /
2) (Γ − / Γ when the cross heat flux isincluded in the model and σ qα = 0 when it is not. Eqns. (18), (19) and (21) may be combined to produceexpressions for ˆ p = ˆ p i + ˆ p e and ˆ p e . Thus ˆ p = − E t ˆ γ ˆ ξ ˆ d i − ˆ c sp ˆ γ i ˆ ∇ · ˆ v + (cid:0) C pe + ˆ c sq (cid:1) ˆ Q − (cid:0) C pe + ˆ c sq (cid:1) i ˆ λ ˆ B x + 2ˆ c sq i ˆ k k ˆ B ′ x , (25)and ˆ p e = − E e ˆ γ ˆ ξ ˆ d i − ˆ c spe ˆ γ i ˆ ∇ · ˆ v + (cid:0) C pe + ˆ c sqe (cid:1) ˆ Q − (cid:0) C pe + ˆ c sqe (cid:1) i ˆ λ ˆ B x + 2ˆ c sqe i ˆ k k ˆ B ′ x , (26)where ˆ c sqe = C qe ˆ c se , ˆ c sq = C qe ˆ c se + C qi ˆ c si , ˆ c sp = ˆ c spe + ˆ c spi , ˆ c spi = ˆ c si ˆ γ i − i ˆ ω ∗ q i / Γˆ γ i − i ˆ ω ∗ qi , (27) ˆ c spe = ˆ c se ˆ γ pe − i ˆ ω ∗ q e / Γˆ γ pe − i ˆ ω ∗ qe , (28) C pe = σ pe i ˆ ω ∗ e − Γ f T e i ˆ ω ∗ n ˆ γ pe − i ˆ ω ∗ qe , (29) E i = (ˆ γ i − i ˆ ω ∗ qi ) − (cid:18) ˆ ω ∗ i − f T i i ˆ ω ∗ n ˆ ω ∗ q i ˆ γ i (cid:19) , (30) E e = (ˆ γ pe − i ˆ ω ∗ qe ) − (cid:18) ˆ ω ∗ e − f T e i ˆ ω ∗ n ˆ ω ∗ q e ˆ γ i − ˆ ω ∗ n f T e Γ σ pe i ˆ ω ∗ ˆ γ i (cid:19) , (31)and E t = E i + E e . IV. SYSTEM OF EQUATIONS
We next algebraically reduce Eqns. (10), (13)-(17), (25) and (26) from a system of eight equations to asystem of five. These five equations use ˆ B x , ˆ ∇ · ˆ v , ˆ Q , ˆ ξ , and ˆ v k as primary variables. Two of these areunmodified from the system of eight: the radial induction equation, Eqn. (10), and the parallel velocityequation, Eqn. (16). One is slightly modified: the parallel induction equation provides an expression for ˆ Q after ˆ n and ˆ p are eliminated. And two new equations are derived: an expression for ˆ ∇ · ˆ v and a parallelvorticity equation which governs ˆ ξ .Eqns. (15) and (16) are combined to provide an expression for ˆ ∇ · ˆ v , ˆ p = ˆ γ i ˆ ∇ · ˆ v − ˆ γ i ˆ γ ˆ ξ ′ − ˆ Q + i ˆ k k ˆ ω ∗ ˆ d i ˆ B r + i (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ⊥ + i ˆ k k (cid:16) ˆ ∇ · ˆ Π gv (cid:17) k . (32)After multiplying by ˆ γ i and substituting Eqn. (25) for ˆ p , ˆ c sp ˆ ∇ · ˆ v = ˆ γ i ˆ γ ˆ ξ ′ − ˆ γ i E t ˆ γ ˆ ξ ˆ d i + ˆ γ i (cid:0) C pe + ˆ c sq (cid:1) ˆ Q − ˆ γ i ˆ k k ˆ ω ∗ ˆ d i ˆ λ + C pe + ˆ c sq ! i ˆ λ ˆ B x + 2ˆ γ i ˆ c sq i ˆ k k ˆ B ′ x − i ˆ γ i (cid:20)(cid:16) ˆ ∇ · ˆ Π gv (cid:17) ⊥ + ˆ k k (cid:16) ˆ ∇ · ˆ Π gv (cid:17) k (cid:21) . (33)The inertial contributions ( ˆ γ i ˆ ∇ · ˆ v ) are dropped as they are small compared to the ˆ c sp ˆ ∇ · ˆ v term from ˆ p inEqn. (25). Without drift and FLR effects only the first and third terms on the right side contribute to ˆ ∇ · ˆ v .The second term on the right side is a drift-like term from ˜ v · ∇ p and ˜ v · ∇ n and the remaining terms arecontributions from electron advection ( ∼ C pe ), cross heat flux ( ∼ c sq ) and ion gyroviscosity.After eliminating ˆ B k , ˆ n and ˆ p from the parallel induction equation, Eqn. (13), we find (ˆ γ i − i ˆ ω ∗ ) ˆ Q = ( A −
1) ˆ ∇ · ˆ v + i ˆ k k ˆ v k + ˆ k k ˆ d i ˆ B ′′ x + (cid:18) ˆ ω ∗ + i ˆ ω ∗ n Γˆ c s E n (cid:19) ˆ γ ˆ ξ ˆ d i − (cid:20) i ˆ ω ∗ + i ˆ ω ∗ n Γˆ c s (cid:0) C pe + ˆ c sqe (cid:1)(cid:21) ˆ Q + i ˆ ω ∗ n Γˆ c s (cid:0) C pe + ˆ c sqe (cid:1) i ˆ λ ˆ B x + S − g ˆ Q ′′ , (34)where A = i ˆ ω ∗ i ˆ γ i + Γ ˆ c spe ˆ c s i ˆ ω ∗ n ˆ γ i , (35)and E n = E e + ˆ ω ∗ i ˆ γ i . (36)Without drift effects, all contributions from ∇ p e and ∇ n vanish (the latter of these results from the /ne factors in Ohm’s law). In particular, these contributions lead to the A , E n , C pe and ˆ c sq factors in Eqn. (34).The only unused equation from our original system of eight is the radial momentum equation, Eqn. (14).To find an expression for ˆ p ′ we take the derivative of Eqn. (32): ˆ p ′ = − ˆ γ i ˆ γ ˆ ξ ′′ − ˆ ω ∗ n Γˆ c s ˆ γ i ˆ γ ˆ ξ ′ ˆ d i + ˆ ω ∗ ˆ d i ˆ Q − ˆ Q ′ + i ˆ k k ˆ ω ∗ ˆ d i ˆ B ′ x + ˆ ω ∗ ˆ d i + ˆ ω ∗ ˆ k k ˆ λ ˆ d i ! i ˆ λ ˆ B x + i (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ′⊥ + i ˆ k k (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ′k + i ˆ λ (cid:16) ˆ ∇ · ˆ Π gv (cid:17) k . (37)Again, we ignore the inertial term ( ˆ γ i ˆ ∇ · ˆ v ). Substituting into Eqn. (14) and applying the tearing ordering, ˆ γ i ˆ γ ˆ ξ ′′ = 2ˆ k k ˆ λ ˆ Q − i ˆ k k ˆ B ′′ x + ˆ ω ∗ n Γˆ c s ˆ γ i ˆ γ ˆ ξ ′ ˆ d i − i ˆ ω ∗ ˆ d i ˆ λ ˆ B x − (cid:16) ˆ ∇ · ˆ Π (cid:17) r − i (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ′⊥ − i ˆ k k (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ′k − i ˆ λ (cid:16) ˆ ∇ · ˆ Π gv (cid:17) k . (38)Without drift and FLR effects, this equation becomes the standard form of the parallel vorticity equation, ˆ γ i ˆ γ ˆ ξ ′′ ≃ − i ˆ k k ˆ B ′′ x .We now have a system of five equations: Eqns. (10), (16), (33), (34), and (38). The discussion of thetearing-ordered contributions from ion gyroviscosity are deferred until the next section. Without thesecontributions, compressibility and parallel flows only couple to this system through the parallel inductionequation, Eqn. (34). Thus in the single-fluid regime where the Hall effect and ion gyroviscosity may beignored, only two equations, the radial induction and parallel vorticity equations, are required for a solution. A. Considerations of Ion Gyroviscosity
With tearing-ordered gyroviscous contributions, the compressibility equation (Eqn. (33)) becomes ˆ c sp ˆ ∇ · ˆ v = ˆ γ i E t ˆ γ ˆ ξ ˆ d i + ˆ γ i (cid:0) C pe + ˆ c sq (cid:1) ˆ Q − ˆ γ i ˆ k k ˆ ω ∗ ˆ d i ˆ λ + C pe + ˆ c sq ! i ˆ λ ˆ B x + 2ˆ γ i ˆ c sq i ˆ k k ˆ B ′ x + ˆ γ i ˆ γ gvi ˆ γ ˆ ξ ′ − σ gv i ˆ γ i c si Γ ˆ d i ˆ γ ˆ ξ ′′ , (39)and the parallel-momentum equation (Eqn. (16)) becomes ˆ γ gvi ˆ v k = − ˆ ω ∗ ˆ d i ˆ B x + i ˆ k k ˆ c sp ˆ γ i ˆ ∇ · ˆ v − i ˆ k k ˆ γ i (cid:0) C pe + ˆ c sq (cid:1) ˆ Q + i ˆ k k E t ˆ γ ˆ ξ ˆ d i − σ gv ˆ c si Γ ˆ λ ˆ d i ˆ γ ˆ ξ ′ − σ gv ˆ c si Γ ˆ d i ˆ k k ˆ γ ˆ ξ ′′ , (40)0where σ gv is a marker for ion gyroviscosity (set to unity when gyroviscosity is included and otherwise zero),the modified ion gyroviscous frequency is ˆ γ gvi = ˆ γ ExB + i ˆ ω ∗ i − σ gv (cid:18) i ˆ ω ∗ i − i ˆ ω ∗ ˆ c si Γ (cid:19) , (41)and ˆ γ ExB is the doppler-shifted growth rate. The tearing-ordered ion-gyroviscous contributions to parallel-vorticity equation (Eqn. (38)) are − (cid:16) ˆ ∇ · ˆ Π (cid:17) r − i (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ′⊥ − i ˆ k k (cid:16) ˆ ∇ · ˆ Π gv (cid:17) ′k − i ˆ λ (cid:16) ˆ ∇ · ˆ Π gv (cid:17) k = i ˆ ω ∗ ˆ ω ∗ i ˆ d i ˆ ∇ · ˆ v − i ˆ ω ∗ (cid:18) ˆ ω ∗ i + ˆ ω ∗ ˆ c si Γ (cid:19) ˆ γ ˆ ξ ′ ˆ d i − i ˆ c si Γ ˆ d i (cid:18)(cid:16) ˆ ∇ · ˆ v (cid:17) ′′ + i ˆ k k ˆ v ′′k (cid:19) − (cid:18) i ˆ ω ∗ i + i ˆ ω ∗ ˆ c si Γ (cid:19) ˆ γ ˆ ξ ′′ . (42)The i ˆ ω ∗ i ˆ γ ˆ ξ ′′ term produces the standard gyroviscous cancellation and cancels the advective diamagneticdrift, however, as there are many additional terms in the this equation, this cancellation is inexact. The i ˆ ω ∗ ˆ c si ˆ γ ˆ ξ ′′ / Γ term is the result of a drift proportional to the gradient of the magnetic field as previouslydiscussed in detail for tearing in a cylindrical pinch configuration [7] (it has been re-characterized in termsof ω ∗ through equilibrium force balance). Combining Eqns. (38) and (42) and again applying the tearingordering gives − ˆ γ gvi ˆ γ ˆ ξ ′′ = 2ˆ k k ˆ λ ˆ Q − i ˆ k k ˆ B ′′ x − i ˆ ω ∗ ˆ d i ˆ λ ˆ B x − iσ gv ˆ c si Γ ˆ d i (cid:18)(cid:16) ˆ ∇ · ˆ v (cid:17) ′′ + i ˆ k k ˆ v ′′k (cid:19) . (43)The last two terms on the right side of Eqn. (43) raise the differential order of the system of equations.Without these contributions, compressibility and parallel flows can be eliminated algebraically from theparallel induction equation, Eqn. (34), which is the only other location where these variables enter thesystem of equations. We do not presently have a solution to the system of equations with ion gyroviscosity,and thus we proceed without the full contributions.Prior work typically includes only the standard gyroviscous cancellation as a model of ion gyroviscosity.Although we can not justify this approximation from a tearing-ordered-equations stand point, we retain the ˆ γ gvi terms as is in order to facilitate comparison. The two relevant limits are then without gyroviscosity( ˆ γ gvi → ˆ γ i ), and with the exact gyroviscous cancellation ( ˆ γ gvi → ˆ γ ExB ). B. Tearing Normalized System of Equations
Without ion gyroviscosity, compressibility and parallel flow can be eliminated algebraically. SubstitutingEqns. (16) and (33) into Eqn. (34) we find ˆ τ Q ˆ Q = ˆ k k ˆ d i ˆ B ′′ x + S − g ˆ Q ′′ − ˆ k k ˆ γ gvi ˆ Q + ˆ τ B − ˆ k k ˆ ω ∗ ˆ γ gvi ˆ d i ˆ λ ! i ˆ λ ˆ B x + ˆ τ ξ ˆ γ ˆ ξ ˆ d i , (44)where ˆ τ Q = ˆ γ i + i ˆ ω ∗ n Γˆ c s (cid:0) C pe + ˆ c sqe (cid:1) − ˆ γ i ˆ c sp (cid:0) C pe + ˆ c sq (cid:1) ( A − , (45) ˆ τ B = i ˆ ω ∗ n Γˆ c s (cid:0) C pe + ˆ c sqe (cid:1) + ˆ γ i (cid:0) C pe + ˆ c sq (cid:1) ( A − c sp , (46)and ˆ τ ξ = ˆ ω ∗ + i ˆ ω ∗ n Γˆ c s E n − ( A − c sp ˆ γ i E t . (47)1Equations (10), (43) (with σ gv = 0 ), and (44) now comprise our system of equations for ˆ B x , ˆ Q and ˆ ξ .The first two terms on the right side of Eqn. (44) are the contributions from the Hall term and resistivediffusion, respectively; the remaining terms result from a combination of compressibility, parallel flows, ∇ p e contributions, inertia and the v × B term. Compressibility and parallel flows contribute the ˆ k k and ˆ ω ∗ termson the right side of Eqn. (44) as well as the ˆ γ i / ˆ c sp terms in the ¯ τ factors. The ∇ p e term in Ohm’s lawcontributes the ˆ ω ∗ n / ˆ c s terms in the ¯ τ factors.With the constant- ψ approximation, where ˆ B x is assumed constant within the small tearing layer, Eqn. (43)is used to eliminate ˆ B ′′ x ; which results in a system of two coupled equations for ˆ Q and ˆ ξ . We use a tearingnormalization for these equations similar to Ref. [8] with the dimensionless variables, ¯ x = ˆ x ˆ d , ¯ ξ = i ˆ k ′k ˆ d ˆ γ ˆ ξ ˆ B r (0) ˆ γ e , ¯ Q = ˆ k ′k ˆ d ˆ d i ˆ Q ˆ B r (0) ˆ γ e , (48)and the dimensionless parameters, ˆ d = ˆ γ ExB ˆ k ′k S g / , ¯ σ = ˆ γ ExB ˆ d i ˆ k ′k ˆ d = ˆ γ ExB ˆ d i S g , ¯ R = ˆ γ gvi ˆ γ ExB , ¯Λ = i ˆ ω ∗ ˆ γ e ˆ γ ExB ˆ γ gvi , (49) ¯ τ Q = ˆ γ ExB (cid:16) ˆ k ′k ˆ d (cid:17) ˆ τ Q , ¯ τ ξ = i ˆ γ ExB (cid:16) ˆ k ′k ˆ d (cid:17) ˆ τ ξ , and ¯ τ B = i ˆ γ ExB ˆ d i ˆ γ e ˆ d (ˆ τ B + 2 i ˆ ω ∗ ) . (50)With this normalization, ¯ σ is the ion skin depth, d i , normalized to the tearing skin depth, δ = ( S g ˆ γ ) − / .Validity of a first-order FLR model requires ρ i < δ . A good rule of thumb for plasmas with comparable ionand electron temperatures is to use the ion sound gyroradius, ρ s = c s /ω ci , and require ρ s /δ = ˆ c s ¯ σ = √ β ¯ σ < .After expanding k k and retaining only the leading order term in x , k ′k x , ¯ R ∂ ¯ ξ∂ ¯ x = ¯ x (cid:0) ¯ ξ + ¯ Q (cid:1) − ¯ x , (51)and ∂ ¯ Q∂ ¯ x = (cid:0) ¯ R − ¯ x + ¯ τ Q (cid:1) ¯ Q + ¯ R ¯ σ ∂ ¯ ξ∂ ¯ x + ¯ τ ξ ¯ ξ − ¯ τ B + ¯Λ¯ x (52)compose the system of second-order coupled equations.Equation (51), a combination of the radial-induction and parallel-vorticity equations, governs the iondynamics and is composed of the contribution from resistivity on the left side, and the contributions fromthe v × B , Hall , and inertial terms, respectively, on the right side. In the single-fluid limit where the Hallterm ( ¯ x ¯ Q ) can be ignored, this equation alone governs the bulk-flow-mediated mode dynamics. Equation(52), a combination of the parallel-induction and parallel-vorticity equations, governs the electron dynamics.The left side of this equation is the contribution from the diffusion of the parallel field and third term on theright side is the contribution from the Hall term ( ˆ k k ˆ d i ˆ B ′′ x ). The ¯ τ parameters scale as β − and are typicallyimportant only at small values of β . The dominant β − contributions result from the gradient of the electronpressure in Ohm’s law (terms involving ˆ ω ∗ n ) and perpendicular compressibility (otherwise). There are othercontributions to ¯ τ Q and ¯ τ ξ from the parallel-field inertia and the v × B term, respectively, however these termare unimportant from a practical perspective. The first and last term on the right side are also contributionsfrom perpendicular compressibility and are important in the moderate- β transition regime (PR2). V. DRIFT-TEARING DISPERSION RELATIONS BY PARAMETRIC REGIME
Once solutions for ¯ Q and ¯ ξ are found, the dispersion relation may be computed by integrating the radialinduction equation (Eqn. (10)) and applying the boundary condition ˜ B ′ r ( ±∞ ) = 0 . The resulting equation2 PR1a PR2 PR3 PR4 PR5 PR1bRegime ¯ σ << σ ∼ σ >> σ >> σ >> τ Q >> ¯ σ boundary and and and and and and ¯ τ Q << τ Q << τ Q << ¯ σ ¯ τ Q ∼ ¯ σ ¯ σ << ¯ τ Q << ¯ σ ¯ τ Q >> or ¯ τ Q << ¯ σ Dominant field ¯ ξ ¯ ξ and ¯ Q ¯ Q ¯ Q ¯ Q ¯ ξ and ¯ QB x diffusion " " " " " " B k diffusion " " " Hall decoupling " " " " ∇ · v decoupling " " " no drift reference [1] [8] [10] [11] [9] [1]drift reference new new new new [9] [6]Hall drift " " " " " " ∇ p e drift " " " ∇ · v drift " " " " Table II: A summary of the parametric regime boundaries, significant terms and fields, and prior references if appli-cable. is D = ˆ ∞−∞ d ¯ x (cid:0) − ¯ x ¯ ξ − ¯ x ¯ Q (cid:1) = ˆ k ′k / ˆ∆ ′ ˆ γ e ˆ γ / ExB S / g (53)where we have defined D for notational convenience. The right side of this expression is the contributionfrom resistivity, thus the integrand of left side of this expression is the ideal radial Ohm’s law. As resistivityis only significant in the layer, proper matching of the inner and outer region solutions ensures the integrandvanishes outside the layer and the integral converges.We next derive the dispersion relation in the various parametric regimes as summarized in Tab. II. Webegin in the single-fluid regime (PR1) with ¯ τ Q << (near PR2) and work our way clockwise around Fig. 1.We do not address PR6, which was solved numerically in Ref. [8], as it is of limited relevance to fusion-plasmaexperiments. We finish again in the single-fluid regime (PR1) with ¯ τ Q >> ¯ σ (near PR6) where we recoverthe drift-tearing result of Ref. [6]. A. PR1a
We use PR1a as a notation for the upper left quadrant of Fig. 1 where ¯ τ Q , ¯ τ ξ , ¯ τ B , ¯ σ , ¯Λ << ∼ ¯ x ∼ ¯ ξ . (54)Examination of the system of tearing equations (Eqns. (51) and (52)) shows ¯ Q << ¯ ξ . Thus the electronequation (Eqn. (52)) may be ignored and the governing equation is simply ¯ R ¯ ξ ′′ = ¯ x ¯ ξ − ¯ x . (55)The solution for ¯ ξ can be expressed in terms of the parabolic cylinder function, U (0 , ¯ x ) = ¯ x ˆ dµ (cid:0) − µ (cid:1) − / exp (cid:20) − µ ¯ x (cid:21) , (56)as ¯ ξ = ¯ R − / U (cid:0) , ¯ R − / ¯ x (cid:1) . Integrating Eqn. (53), the drift dispersion relation is ˆ γ e ˆ γ / gvi = ˆ γ / MHD (57)3where ˆ γ MHD is the single-fluid growth rate without drift effects, ˆ γ MHD = S − / g ˆ∆ ′ √
2Γ ( / ) ! / ˆ k ′k / . (58) B. PR2
Regime PR2 is the transition at moderate β between the single-fluid regime, PR1, and the electron-MHDregime, PR3. Here we assume ¯Λ ∼ ¯ x ∼ , ¯ ξ ∼ ¯ Q and ¯ τ Q , ¯ τ ξ , ¯ τ B << or ¯ τ Q , ¯ τ ξ , ¯ τ B << ¯ σ . (59)Thus the system of tearing equations becomes ¯ Q ′′ = ¯ R − ¯ x ¯ Q + ¯ R ¯ σ ¯ ξ ′′ + ¯Λ¯ x , (60)and ¯ R ¯ ξ ′′ = ¯ x (cid:0) ¯ Q + ¯ ξ (cid:1) − ¯ x . (61)Following the method outlined in Ref. [8] for the solution of a similar system of equations (where ¯ R → and ¯Λ → ), we transform this system of equations into two independent parabolic cylinder equations, λ − i ¯ V ′′ i = ¯ x ¯ V i − C i ¯ x , (62)where ¯ V i = ¯ ξ + a i ¯ Q and i = 1 , . This transformation requires λ i = ¯ R − + ¯ σ a i , (63) a i = 12 ± r R ¯ σ , (64)and C i = 1 − ¯Λ ¯ R ( a i − . (65)The solution for each ¯ V i is ¯ V i = λ / i C i U (cid:16) , λ / i ¯ x (cid:17) . Integrating Eqn. (53) to find the dispersion relationgives D = 2 π Γ ( / )Γ ( / ) " C a λ − / − C a λ − / a − a . (66)This may be expressed in a more explicit form as D = √
2Γ ( / ) f (cid:0) ¯ σ, ¯ R, ¯Λ (cid:1) , where f (cid:0) ¯ σ, ¯ R, ¯Λ (cid:1) = X i =1 , " − ¯Λ ¯ R ( − i r R ¯ σ − ! − i (cid:18) R ¯ σ (cid:19) − / × " ¯ R − + ¯ σ − i ¯ σ r ¯ σ R − − / . (67)The limits of this expression under the same approximations as PR1a and PR3 are consistent with thedispersion relations found in these regimes. Consider the limit where ¯ σ << , in this case f (cid:0) ¯ σ, ¯ R, ¯Λ (cid:1) → (cid:0) R/ (cid:1) ¯ R / . With the additional limit ¯Λ << (as is the case in PR1a), f (cid:0) ¯ σ, ¯ R, ¯Λ (cid:1) → ¯ R / and werecover Eqn. (57). In the limit where ¯ σ >> , f (cid:0) ¯ σ, ¯ R, ¯Λ (cid:1) → ¯ σ − / . As we shall see in the next subsection,this limit is the dispersion relation found in the electron-MHD regime, PR3.4 C. PR3
In the electron-MHD regime, the resistive diffusion of B k balances the Hall term in the parallel inductionequation, and the parallel-vorticity equation is not needed. The orderings of this regime are a small tearinglayer and large B k , ¯ x − ∼ ¯ σ / ∼ ¯ Q , small ion displacement, ¯ ξ ∼ ¯ σ − / , and large d i , ¯ σ >> , such that ¯Λ << ¯ σ , ¯ τ ξ << ¯ σ , ¯ τ B << ¯ σ / and ¯ τ Q << ¯ σ . After substituting the ordered electron equation, Eqn. (52),into the ordered ion equation, Eqn. (51), the governing equation in this regime is ¯ Q ′′ = ¯ σ ¯ x ¯ Q − ¯ σ ¯ x . (68)The solution to this equation is ¯ Q = √ ¯ σU (cid:0) , √ ¯ σ ¯ x (cid:1) . Integrating Eqn. (53), we find D = √
2Γ ( / ) ¯ σ − / (the limit of D from PR2 when ¯ σ >> ) and the dispersion relation is then ˆ γ e = S − / g (cid:16) ˆ d i ˆ k ′k (cid:17) / ˆ∆ ′ √
2Γ ( / ) . (69)In this regime the growth rate scales as d / i S − / and the mode simply rotates at the electron drift frequency;there is no drift stabilization. This result is not particularly surprising, as the mode is mediated purely bythe electron fluid through the induction equation. Contributions from ion compressibility, parallel ion flowsand ion vorticity do not play a role. D. PR4
The PR4 regime is the transition between the B k -diffusion (PR3) and the semicollisional (PR5) regimes.The orderings of this regime are similar to PR3; a small tearing layer with a large B k , ¯ x − ∼ ¯ σ / ∼ ¯ Q ,small ion displacement, ¯ ξ ∼ ¯ σ − / , however ¯ τ Q is comparable to the normalized ion skin depth which islarge, ¯ σ >> , such that ¯Λ << ¯ σ , ¯ τ ξ << ¯ σ , ¯ τ B ∼ ¯ σ / and ¯ τ Q ∼ ¯ σ . Thus the ¯ τ Q and ¯ τ B contributionsmust both be retained in Eqn. (52), and the system of tearing equations becomes ¯ Q ′′ = ¯ τ Q ¯ Q + ¯ R ¯ σ ¯ ξ ′′ − ¯ τ B , (70)and ¯ R ¯ σ ¯ ξ ′′ = ¯ σ ¯ x ¯ Q − ¯ σ ¯ x . (71)These equations may be combined into a single non-homogeneous parabolic cylinder equation for ¯ Q , ¯ Q ′′ = ¯ τ Q ¯ Q + ¯ σ ¯ x ¯ Q − ¯ σ ¯ x − ¯ τ B . (72)The solution to this equation up to a constant of integration, following the method outlined in Ref. [18], is ¯ Q (¯ x ) = A + ˆ ∞ exp ik ¯ x r ¯ σ ! U ( a, k ) U ( a, dk + A − ˆ − ∞ exp ik ¯ x r ¯ σ ! U ( a, k ) U ( a, dk , (73)with the constraint A + − A − = i ¯ σ / √ τ B σ U ( a, U ′ ( a, , (74)where a = ¯ τ Q / σ . The constant of integration (either A + or A − ) is found by matching the layer equationswith the outer solution (in practice, requiring that the integral of Eqn. (53) converges), which provides theadditional condition A + = − i p ¯ σ/ . Integrating Eqn. (53) determines the dispersion relation as ˆ γ e Γ [(3 + ¯ τ Q / ¯ σ ) / τ Q / ¯ σ ) /
4] = S − / g q ˆ k ′k ˆ d i ˆ∆ ′ π . (75)5Drift effects modify the dispersion relation through the left side, in particular the drift modified growth rate ˆ γ e and the drift effects contained in ¯ τ Q and ¯ σ . In the limit where ¯ τ Q << ¯ σ , the left side of Eqn. (75) becomes ˆ γ e Γ ( / ) / √ π , consistent with the dispersion relation of PR3. In the opposite limit, where ¯ τ Q >> ¯ σ the leftside of the equation becomes ˆ γ e √ ¯ τ Q / √ ¯ σ , which is consistent with the dispersion relation found in the nextsection for the semicollisional regime, PR5. Although ¯ τ B , which scales similarly in magnitude to ¯ τ Q , affectsthe eigenfunction, it does not modify the growth rate. In the limit of PR3, both ¯ τ B and ¯ τ Q are small andthus the results are consistent. In the limit of PR5, where ¯ τ B is again expected to be large, ¯ τ B contributesan even parity term to the eigenfunction and thus again does not contribute to the dispersion relation afterintegration of Eqn. (53). E. PR5
The orderings in the semicollisional regime are similar to PR3 and PR4, with a large B k , ¯ x − ∼ ¯ σ / ∼ ¯ Q ,and small ion displacement, ¯ ξ ∼ ¯ σ − / . However in this regime the ¯ τ terms are larger than normalizedion skin depth (but not too large), ¯ σ >> , such that ¯Λ << ¯ σ , ¯ τ ξ << ¯ σ , and ¯ σ << ¯ τ Q ∼ ¯ τ B << ¯ σ .The diffusion of B k may be neglected and the Hall term in Eqn. (52) is balanced by the ¯ τ Q and ¯ τ B terms.After substitution of the ordered electron equation, Eqn. (52), into the ordered ion equation, Eqn. (51), thegoverning equation for this regime is τ Q ¯ Q + ¯ σ ¯ x ¯ Q − ¯ σ ¯ x − ¯ τ B . (76)The solution is algebraic, ¯ Q = ¯ σ ¯ x + ¯ τ B ¯ σ ¯ x + ¯ τ Q . (77)The dispersion relation, found by integrating Eqn. (53), is then ˆ γ e ˆ τ / Q = S − / g ˆ k ′k ˆ∆ ′ π ˆ d i . (78)The growth rate scales as ρ / s S − / g and drift effects are contained on the left side of Eqn. (78). Theeigenfunction, ¯ Q , only contains even-parity contributions from ¯ τ B , which vanish during the integration ofEqn. (53) and thus do not contribute to the dispersion relation.This is the two-fluid drift-regime first described by Drake and Lee [9]. Simplifying this expression furtherby assuming β << (which defines this regime), σ pe = σ qi = − σ qe = 1 , we find ˆ γ / e (cid:18) ˆ γ ExB ˆ γ i f T i ˆ γ e + f T e ˆ γ i (cid:19) / = S − / g ˆ k ′k ˆ∆ ′ π ˆ c s ˆ d i ! / . (79)When the electron temperature is much larger than the ion, f T i = 0 and f T e = 1 , the standard two-thirds,one-third dispersion relation is attained. Our results are not identical to Drake and Lee; however we donot include ion and electron gyroviscosity or heat-flux contributions to the frictional force. The inclusion ofcross heat flux cancels contributions from the pure density-gradient drifts in the dispersion relation, as with σ qα = 0 but σ pe = 1 , one instead finds ˆ γ / e (cid:18) (ˆ γ ExB − f T e Γ i ˆ ω ∗ n ) (ˆ γ ExB + f T i Γ i ˆ ω ∗ n )ˆ γ e (cid:19) / = S − / g ˆ k ′k ˆ∆ ′ π ˆ c s ˆ d i ! / . (80) F. PR1b
This regime is the low- β , drift limit of the single-fluid regime. With the corresponding orderings, << ¯ τ Q ∼ ¯ τ ξ , ¯ σ << ¯ τ Q ∼ ¯ τ ξ and ¯ τ B ∼ ¯Λ ∼ , the ¯ τ Q and ¯ τ ξ terms balance in the electron equation, Eqn. (52),6and thus ¯ ξ ∼ ¯ Q as ¯ τ Q ¯ Q = − ¯ τ ξ ¯ ξ . Substituting this balance into Eqn. (51) the governing equation becomes ¯ R ¯ ξ ′′ = ¯ x (cid:18) − ¯ τ ξ ¯ τ Q (cid:19) ¯ ξ − ¯ x . (81)The solution is ¯ ξ = ¯ R − / M − / U (cid:0) , ¯ R − / M / ¯ x (cid:1) where M = 1 − ¯ τ ξ / ¯ τ Q . Assuming β << , and thus M ≃ ˆ γ e / ˆ γ i , integration of Eqn. (53) gives the dispersion relation, ˆ γ / MHD = ˆ γ / e ˆ γ / i ˆ γ / gvi . (82)With an exact gyroviscous cancellation, ˆ γ gvi = ˆ γ ExB , this is the standard drift-tearing dispersion relation asfound by Coppi [6].
VI. VERIFICATION OF THE NIMROD CODE
We may now use the dispersion relations of Sec. V to verify the implementation of the unreduced extended-MHD equations in the initial-value NIMROD code [13]. NIMROD is primarily designed as a nonlinear-physicscode. However, it uses the linear response of the perturbed system as a preconditioner during nonlinear solves.This functionality makes available the linearized equation within NIMROD and thus permits our verificationexercise. This verification is a partial test of the NIMROD equation implementation as well as a test of thetime and spatial discretizations.Cases are implemented as a periodic-in-y, symmetric-in-z box within NIMROD. Each specific equilibriumis generated by specifying the equilibrium magnetic-shear-scale length ( L S ), the ratio of the magnetic shearto guide field ( ǫ B ), the plasma β , the equilibrium pressure-gradient-scale length ( L P ), and the ratio of thesheared to background pressure ( ǫ P ). Equilibrium fields are computed by solving the MHD-force balanceequations based on a hyperbolic-secant-squared parallel-current profile, λ = µ J · B B = ǫ B L S sech (cid:18) xL S (cid:19) , (83)and a hyperbolic-tangent pressure profile, Γ µ p B = β (cid:18) ǫ P tanh (cid:18) xL P (cid:19)(cid:19) . (84)Our cases use comparable magnetic-shear and pressure-gradient scale lengths, L s = L p , and impose thisgradient with a dominant density profile to avoid ITG-like modes (see Ref. [19]). The fraction of the pressuregradient that results from the density profile, f n = n ′ p /n p ′ , always equals or exceeds / . Drift effectsare included when ǫ P = 0 . For cases with ǫ P = 0 , the tearing stability parameter, ∆ ′ , may be computedanalytically for this equilibrium (Ref. [8]): ∆ ′ = 2 L S (cid:18) kL S − kL S (cid:19) . (85)For cases with ǫ P = 0 , we use NIMROD to infer that ∆ ′ is unchanged. As p ′ is increased, if the growthrate from NIMROD computations with a single-fluid model is unchanged then ∆ ′ is constant. Equation(85) assumes an infinite-in-x domain. This is, of course, not practical for the NIMROD finite-elementcomputations where instead a large ratio of D x /L s is used to approximate the infinite domain, where D x isthe box half length in the x dimension. Our cases use D x /L s = 6 with a 96 radial bi-cubic elements packednear the resonant surface where the single-fluid growth rate discrepancy between NIMROD and the analyticsis less than .Table III summarizes the parameters used for our verification studies. In a practical verification exercise,the physical parameter space (equilibrium characteristic values, length scales and gradients) affect the derivedparameter space (Lundquist number, tearing stability parameter, ion skin depth, β and drift frequencies)in a complex manner. The locations of these cases in the ¯ σ − ¯ τ Q parameter space in the limit where ˆ ω ∗ → is superimposed onto Fig. 1. In general as ˆ ω ∗ increases, ¯ σ marginally increases and the ¯ τ parameters7 case k ⊥ d i β ¯ σ ¯ τ Q ¯ τ B ¯ τ ξ ¯Λ regime stabilizationA (Fig. 2) .
002 0 . . × − .
064 0 . .
048 0 .
65 0.37 0.061 0.29 066 PR3/PR4 noneD (Fig. 6) .
048 1 . × −
54 22 0.017 18 0.87 PR4-PR6 moderateTable III: Parameters for the verification ω ∗ scans with the NIMROD code. The normalized parameters are evaluatedat ˆ ω ∗ = 1 . × − and are modified by drift effects. All scans use σ pe = σ qi = − σ qe = 1 , S = 3 . × , ˆ∆ ′ = 1 . , kL s = 0 . and ǫ B = 0 . . Cases use the electron-to-ion-mass ratio from a Deuterium gas discharge of . × − unless otherwise mentioned. PR2 PR0
Figure 2: Scan A growth rates and normalized parameters with v E × B = − v ∇ p and f n = 1 . The converged resultsfrom NIMROD runs (points) are compared with the drift analytics of PR2 (lines, Eqns. (53) and (66)). increase linearly moving the cases down (and slightly to the right) in the ¯ σ − ¯ τ Q parameter space of Fig. 1,as illustrated with dashed lines. Our choice of scan locations in the ¯ τ Q - ¯ σ phase space is the result of acombination of finding a representative sample of cases to fill the experimentally relevant parameter space ofTable I, choosing cases which are able to achieve reasonable ˆ ω ∗ ∝ ǫ P β ˆ d i with ǫ P < (which avoids negativepressure regions), and testing the analytics in a variety of regimes.All cases rotate in the electron diamagnetic direction. The dominant ω ∗ e influence results from thedenominator of the right side of Eqn. (53). In the electron-MHD regime of PR3, where the ion dynamics nolonger influence the mode, the mode is at rest in the frame of the electron fluid.The scan A growth-rate comparison at moderate β (0.1) and low ˆ d i (0.002) between NIMROD runs andthe dispersion relation of PR2 is shown in Fig. 2. Good agreement is achieved until ¯ τ Q ∼ (the five left-mostpoints in the figure agree with the analytics with less than a error) and the mode enters the regime ofPR0. Although there are no analytics for this regime, we note NIMROD predicts stronger drift-stabilizationin PR0 than the relatively weak effect predicted by the drift analytics in PR2. In fact, at larger values of ˆ ω ∗ NIMROD predicts complete stabilization in PR0 as NIMROD cases at ˆ ω ∗ = 7 . × − are stable.Figure 3 shows the scan B result of a verification scan at moderate ˆ d i (0.064) and β (0.1) which again beginsin PR2 and transitions to PR0. Similar to scan A, as ˆ ω ∗ is increased ¯ τ Q approaches unity and the mode8 PR2 PR0
Figure 3: Scan B growth rates and normalized parameters with v E × B = 0 and f n = 0 . . The converged results fromNIMROD runs (points) are compared with the drift analytics of PR2 (lines, Eqns. (53) and (66)). ultimately enters PR0 where there is no analytic solution. However, unlike scan A, both the computationsand the analytics predict complete stabilization of the mode at qualitatively similar values of ω ∗ (NIMRODcomputations at ˆ ω ∗ = 5 . × − are stable). Similar to scan A, the first five left-most points are within 3%of the analytic results.Figures 4 and 5 show the scan C growth-rate comparisons at large ˆ d i (2.048) and moderate β (0.1). Thesescans begin in PR3 and transition into PR4. The NIMROD cases agree within 5% and 1% of the analyticresults for Figs. 4 and 5, respectively. The cases in Fig. 4 are essentially collisionless and there is no driftstabilization as ˆ ω ∗ increases, instead the mode growth rate increases. In the collisionless regime withoutthe advective term in electron inertia, as currently implemented in NIMROD, S g → ˆ γ − ˆ d − e and Eqn. (69)becomes ˆ γ e ˆ γ − = ˆ d e ˆ d i ˆ k ′k ˆ∆ ′
2Γ ( / ) ≡ ˆ γ c . (86)In the limit of this equation where ˆ ω ∗ e << ˆ γ c , the mode grows at the drift-free growth rate and drifts at theelectron drift frequency, ˆ γ ≃ ˆ γ c + i ˆ ω ∗ e . In the limit where ˆ ω ∗ e >> ˆ γ c , the mode grows proportionally to thesquare root of the drift frequency ˆ γ ≃ (cid:0) ˆ γ c + √ ˆ ω ∗ e ˆ γ c (cid:1) / i (cid:0) ˆ ω ∗ e + √ ˆ ω ∗ e ˆ γ c / (cid:1) . This second limit explains thedestabilization of the mode as seen in the figure. It is of interest to note that if the advective term is includedin electron inertia then S g → ˆ γ − e ˆ d − e and there is no growth rate increase. However, electron gyroviscosityenters the equations on the same order and should also be retained. The relevant physical effects within theboundary layer for these near-collisionless cases illustrate the breakdown of the argument to ignore electronadvection and gyroviscosity presented in Sec. II. As the resonant condition causes the dominant terms inOhm’s law to vanish, the boundary layer physics is determined by a balance of the remaining, otherwise smallterms. For the collisionless-drift-tearing mode, these small terms include electron advection and gyroviscosity(Ref. [15] includes these terms in PR5 without drift effects). The cases in Fig. 5 are identical to those in Fig. 4except they are collisional through the use of a small electron mass, µ = 2 . × − . For these collisionalcases the mode is not drift stabilized and simply rotates with the electron fluid as predicted by Eqn. (69)when S g → v A µ /k ⊥ η .9 PR3 PR4
Figure 4: Scan C with v E × B = − v ∇ p and f n = 1 . The converged results from NIMROD runs (points) are comparedwith the drift analytics from PR4 (lines, Eqn. (75)). PR3 PR4
Figure 5: Scan C with a reduced electron mass, µ = 2 . × − , v E × B = 0 and f n = 1 . The converged results fromNIMROD runs (points) are compared with the drift analytics from PR4 (lines, Eqn. (75)). PR4
PR5 PR6
Figure 6: Scan D with µ = 2 . × − , v E × B = − v ∇ p and f n = 0 . . The converged results from NIMROD runs(points) are compared with the drift analytics from PR4 (lines, Eqn. (75)). Beyond the successful verification of the code in this electron-fluid-mediated regime, the validity of themodel remains in question. For first-order electron-FLR model validity, one requires that ρ e /δ = √ βµ ¯ σ << ;a condition that is satisfied for these cases. However, it is unlikely that the simple electron-response modelis sufficient to model the collisionless dynamics of Fig. 4. Given that the ion gyroviscous cancellation isincomplete (see Sec. IV A), the implicit assumption in the model that ∇ · Π e,gv + m e v e · ∇ v e = 0 is likely notvalid when ˆ ω ∗ is large. Further study and code development pertaining to this issue is required and outsidethe scope of this work.The scan D verification exercise that begins in PR4 and transitions through PR5 to PR6 at low β ( . × − ) and large ˆ d i (2.048) is shown in Fig. 6. Although we are not able to run a drift-verification scanwhile starting in the semicollisional regime, PR5, this comparison does include cases near this regime. Inthis regime the mode is weakly stabilized where the growth rate is decreased by approximately a factor offive for large values of ω ∗ /γ . The discrepancy between the analytics and the numerics for the first six casesis approximately 15%, however, the right-most three cases, where the drift effects are large, agree with theanalytic theory within 7%, 2% and 0.2%, respectively.Figure 7 is a matrix of eigenfunction plots for scans A-D at small and large values of ˆ ω ∗ . The scalings ofthe ¯ ξ and ¯ Q are consistent with the assumptions for the various regimes made in Sec. V. For the small- ω ∗ ,scan-A case ¯ ξ >> ¯ Q which is reasonable for a case near the single-fluid limit. When ω ∗ is large (scan A andall of scan B), ¯ ξ ∼ ¯ Q in line with the assumptions of PR2. For cases in PR4 through PR6 (scans C andD), ¯ Q >> ¯ ξ , ¯ Q is larger than unity, and the eigenfunction is more localized consistent with the orderingsof ¯ Q ∼ ¯ σ / ∼ ¯ x − . All cases except the large- ω ∗ , scan-D case produce an odd eigenfunction (only theodd component contributes to the growth rate, a result of Eqn. (53)). The large- ω ∗ scan D case has aneven component which is in agreement with the discussion of Secs. V D and V E and large contributionsfrom ¯ τ B . Finally, at large ω ∗ only scan D has radial drift structures that extend to ¯ x = ± (not shown).All other cases do not exhibit this structure and the eigenfunction is highly localized within the resonantlayer. Unlike previous computational verification drift-tearing work [12], our computations do not exhibitsignificant influence from the computational boundary condition.1 Scan D (cid:0) = 1.3x10 -6 * Scan C (cid:2) = 1.3x10 -6 * Scan A (cid:3) = 1.3x10 -6 * Scan D (cid:4) = 3.3x10 -4 * Scan B (cid:5) = 1.3x10 -6 *^^^^ ^ Scan C (cid:6) = 3.3x10 -4 *^ Scan B (cid:7) =4.5x10 -5 *^ Scan A (cid:8) = 4.5x10 -5 *^ Figure 7: Eigenfunctions ¯ ξ and ¯ Q for scans A-D (top to bottom) at small (left) and large (right) values of ˆ ω ∗ . Theplots for scan C correspond to Fig. 5 VII. CONCLUDING DISCUSSION
This work is both an analytic and computational investigation of drift tearing with an unreduced, extended-MHD model. Our new analytic results have been used to verify the implementation of the extended-MHDequations within the NIMROD code. As the tearing-layer dynamics result from the balance of otherwisesmall terms, this verification is a novel way to test the extended-MHD implementation. Our new analyticresults describe the experimentally relevant portion of the drift-tearing phase space. Within this phasespace, there is the potential for varying degrees of drift stabilization: there is a weakly stabilizing effect ateither small d i and moderate β or at large d i and small β , complete stabilization is possible at moderate2 d i and β and there is no stabilization at large d i and moderate β where the ion dynamics are decoupledfrom the mode. We emphasize that our definition of moderate β encompasses the values that are pertinentfor a fusion reactor ( β ∼ − ). There are some caveats to the applicability of this work when oneconsiders the validity of the first-order ion FLR model. However, we argue that this model may still bequalitatively valid when the ion gyroradius is no longer small, as the mode transitions to one dominatedsolely by the electron-fluid dynamics (given a sufficient electron-dynamics model). Our results can not bedirectly applied tokamak discharges, as we do not retain the effects of ion gyroviscosity and plasma shapingand curvature. Instead, the ultimate benefit of this work is to provide enhanced confidence in nonlinear,extended-MHD, boundary-layer-dynamics computations of tokamak discharges with reconstructed profilesand realistic geometry. Acknowledgments
The authors would like to thank Carl Sovinec and Chris Hegna for stimulating discussions. This workis supported by US Department of Energy grants DE-FC02-06ER54875 and DE-FG02-08ER54972. Thisresearch used resources of the National Energy Research Scientific Computing Center, which is supportedby the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. [1] H. P. Furth, J. Killeen, and M. N. Rosenbluth, Physics of Fluids , 459 (1963).[2] S. I. Braginskii, Transport Properties in a Plasma in Review of Plasma Physics , edited by M. A. Leontovich,Vol. 1 (Consultants Bureau, New York, 1965).[3] P. J. Catto and A. N. Simakov, Physics of Plasmas , 90 (2004).[4] J. J. Ramos, Physics of Plasmas (1994-present) , 082502 (2010).[5] J. J. Ramos, Physics of Plasmas (1994-present) , 102506 (2011).[6] B. Coppi, Physics of Fluids , 1501 (1964).[7] J. R. King, C. R. Sovinec, and V. V. Mirnov, Physics of Plasmas , 042303 (2011).[8] E. Ahedo and J. J. Ramos, Plasma Physics and Controlled Fusion , 055018 (2009).[9] J. F. Drake and Y. C. Lee, Physics of Fluids , 1341 (1977).[10] S. V. Bulanov, F. Pegoraro, and A. S. Sakharov, Physics of Fluids B , 2499 (1992).[11] V. V. Mirnov, C. C. Hegna, and S. C. Prager, Physics of Plasmas , 4468 (2004).[12] D. Biskamp, Nuclear Fusion , 1059 (1978).[13] C. Sovinec, A. Glasser, T. Gianakon, D. Barnes, R. Nebel, S. Kruger, D. Schnack, S. Plimpton, A. Tarditi, andM. Chu, Journal of Computational Physics , 355 (2004).[14] C. Sovinec and J. King, Journal of Computational Physics , 5803 (2010).[15] R. Fitzpatrick, Physics of Plasmas , 042101 (2010).[16] J. J. Ramos, Physics of Plasmas , 112301 (2005).[17] E. Ahedo and J. J. Ramos, Physics of Plasmas (1994-present) , 072519 (2012).[18] R. Hazeltine and J. Meiss, Plasma Confinement , Dover Books on Physics (Dover Publications, 2013).[19] D. D. Schnack, J. Cheng, D. C. Barnes, and S. E. Parker, Physics of Plasmas (1994-present)20