Parasitic momentum flux in the tokamak core

Tokamak plasmas without applied torque routinely rotate spontaneously in the toroidal (symmetry) direction, exhibiting nonzero, sheared toroidal rotation profiles.¹ This so-called “intrinsic” rotation is not only of fundamental interest: toroidal rotation helps suppress certain instabilities² and its shear may reduce turbulent heat transport.³ These advantages are important for future burning plasma devices such as ITER, in which the dominant α-heating will not exert toroidal torque, unlike the neutral beam heating typical of present-day devices.⁴ 

Although experimentally measured intrinsic rotation profiles are very diverse, many exhibit three distinct radial regions: an edge region with co-current rotation (toroidal rotation in the direction of the plasma current $Ip$⁠), a mid-radius “gradient region” where the rotation either becomes increasingly counter-current with decreasing radius (counter-current peaking) or stays relatively flat, and a flat or weakly co-current-peaked central region affected by sawtoothing.^5–9 Previous theoretical,¹⁰ numerical,¹¹ and experimental⁶ work suggests that the edge rotation is driven by the interaction of passing-ion drift orbit excursions with the spatial variation of the turbulent fluctuations. The present work focuses on the “gradient region” at intermediate radius, where the radial variation of plasma parameters is much slower, allowing other effects to compete with those of orbit excursions.

Over the last decade, intrinsic rotation at mid-radius has undergone intense theoretical and experimental investigation. Nonaxisymmetric magnetic fields can strongly affect the toroidal rotation. The present work will focus exclusively on the case of axisymmetric confining magnetic field, for which the conservation of toroidal angular momentum¹² excludes the possibility of a self-generated torque. Intrinsic rotation must therefore result from a nondiffusive component to the momentum flux. Neoclassical (collisional) momentum fluxes are much too small to explain the experimental observations, implying that turbulent transport is dominant.¹ A number of turbulent calculations suggest the presence of a momentum pinch, a component of momentum flux that is proportional to the toroidal rotation itself, rather than its gradient.¹³ However, these models cannot explain the common observation of sheared velocity profiles passing through zero.^5–9 Such measurements imply the presence of a “residual stress,” meaning a momentum flux contribution that is independent of both toroidal rotation and its radial gradient. For up-down symmetric geometries, often a good approximation for tokamak core plasmas, symmetry arguments restrict the leading-order momentum flux terms from driving residual stress.¹⁴ Theoretical work has accordingly focused on symmetry-breaking mechanisms¹⁴ such as $E×B$ shear,¹⁵ up-down-asymmetric geometry,¹⁶ turbulence spreading¹⁷ and intensity gradient,¹⁸ magnetic shear,¹⁹ and the effect of neoclassical flows.²⁰ Particularly challenging to theory are the experimental observations of rotation reversals in the gradient region, in which counter-current rotation peaking suddenly flattens or switches to weak co-current peaking when plasma density or current cross threshold values.^6–9 The rapidity of these reversals suggests that the peaking is not determined by neoclassical flows or other quantities that vary smoothly with plasma parameters, following instead from properties of the turbulence itself, which may suddenly change character, for e.g., as an instability threshold is crossed. Reference 21 quasilinearly evaluated the residual stress due to several higher-order gyrokinetic terms, including a geometrical correction to the $E×B$ drift. In this letter, I identify a correspondence of forms between that residual stress and a free-energy transfer term key to ion Landau damping, which causes free-energy flows within the turbulence to drive a robust, fully nonlinear symmetry-breaking momentum flux. This flux causes counter-current core rotation peaking consistent with experimental profiles, and explains several observations related to rotation reversals.

To develop intuition, consider first a low-frequency axisymmetric density perturbation, as sketched in Fig. 1. At low frequencies and large scales, electron parallel force balance ensures that the nonzonal electrostatic potential $ϕ̃$ is proportional to the nonzonal ion gyrocenter density $ñi$⁠. The pressure gradient and electric field then cause ions to flow out of the dense region along the magnetic field. The poloidal electric field also causes a radial $E×B$ drift that advects counter- (co-)current ion momentum inward (outward), regardless of the signs of $Ip$ and the toroidal magnetic field $BT$⁠.

Key to this mechanism is a dual role for the weak electric field caused by the poloidal variation of the potential $ϕ$ on length scales comparable to the minor radius $r$⁠. The nonvanishing parallel component of this electric field allows it to cause local ion acceleration, resulting in energy transfer between electrostatic potential (⁠$ϕ$⁠) and parallel ion flow (⁠$u∥i$⁠). Because the background plasma gradients predominantly supply energy to even moments of the distribution function (such as density), while odd moments (such as $u∥i$⁠) are subjected to dissipation,^22,23 steady-state energy balance often requires a net transfer of free energy from the potential (a function of even moments) to the ion parallel flows, causing a statistical symmetry breaking in the corresponding energy transfer term. Although toroidal angular momentum conservation does not allow the self-generated electric field to impart a net torque to the plasma, the weak radial $E×B$ drift due to the poloidally varying $ϕ$ may transport toroidal angular momentum in the radial direction. The correlations between the ion parallel flows and the weak radial $E×B$ drift, resulting from the statistical symmetry breaking due to energy transfer, cause this part of the momentum flux to have a preferred sign, independent of the plasma rotation and its radial gradient. In this letter, we will consider this residual stress in two separate cases: first a simpler special case with axisymmetric fluctuations, where the momentum flux occurs due to damping of geodesic acoustic modes (GAMs) via ion parallel flows, and later the general case including nonaxisymmetric fluctuations, where the momentum flux can occur for any turbulent fluctuations in which energy transfer from the potential to ion parallel flow is nonnegligible.

Both calculations use the simplest model capturing the relevant physics: the large-aspect-ratio limit of the electrostatic, isothermal gyrofluid equations in a radially thin geometry,^22,24 written in cgs units as

with species subscript s meaning ions i or electrons e; species charge state $Z$ (−1 for electrons), mass $ms$⁠, and (constant) temperature $Ts0$⁠; fluctuating $ns$ and equilibrium $ns0$ species density (assuming $Zni0=ne0$⁠); $E×B$ drift $uEs≐(c/B)b̂×∇ϕG$⁠; gyroaveraged potential $ϕG≐Γ1sϕ$⁠; curvature operator $K≐−(2c/B)b̂×∇ ln B·∇$ capturing the magnetic drifts and $E×B$ divergence; parallel gradient $∇∥≐b̂·∇$⁠; parallel flow velocity $u∥s$⁠; and dissipation operator $D∥s$⁠. The gyroaveraging operators $Γ0s$ and $Γ1s$ take the low-$k⊥$ limits $(1−Γ0s)→−ρs2∇⊥2$ and $Γ1s→1$⁠, for $ρs=vts/|Ωcs|$ the species gyroradius, with thermal speed $vts≐(Ts0/ms)1/2$ and gyrofrequency $Ωcs≐ZeB0/msc$⁠. We take the safety factor $q=BTr/BpR0$ order unity, so the poloidal field and inverse aspect ratio are comparably small, $Bp/BT∼r/R0≪1$⁠, allowing explicit appearances of $B$ and $R$ to be replaced with representative constants $B0$ and $R0$⁠, and setting $bT≐b̂·ζ̂→±1$ for magnetic and toroidal directions $b̂$ and $ζ̂$⁠. Since the toroidal component of $uEs$ is small in $Bp/BT$⁠, Eqs. (1)–(3) conserve a simplified toroidal angular momentum involving only the zonal (flux-surface) average $〈⋯〉$ of $u∥i$⁠, assuming $〈D∥i〉=0$

with toroidal angular momentum density $Lζ≐mini0bTR0u∥iand flux Πζ≐mini0bTR0(uEix−2Ti0ZeKx)u∥i$⁠, with $uEix≐uEi·∇x, Kx≐K(x)$⁠, and radial (flux-surface) label $x$⁠. Eq. (4) shows that the toroidal angular momentum is advected by the $E×B$ and magnetic drifts, without sources or sinks.

A nondiffusive momentum flux as in Fig. 1 may be driven by geodesic acoustic mode (GAM) damping, which we may treat in a shearless simple-circular geometry, $Kx→(2cbT/B0R0) sin θ$ for poloidal angle $θ$⁠. Following Ref. 25, we retain only one axisymmetric Fourier component from each of Eqs. (1)–(3), specifically $nis≐〈ni sin θ〉, u∥c≐〈u∥i cos θ〉$⁠, and $uEz≐〈uEiθ〉$ with $uEiθ≐uEi·θ̂=bTcB0∂xϕ$⁠. We neglect electron polarization and take low-$k⊥$ gyroaveraging and quasineutrality $Z〈ni sin θ〉≈〈ne sin θ〉$⁠, and electron adiabatic response $Te0〈ne sin θ〉≈ene0〈ϕ sin θ〉$⁠, obtaining

in which $bp≐b̂·θ̂≪1$ and $Ta≐Ti0+ZTe0$⁠. We have taken $〈D∥i cos θ〉→ν∥u∥c$ for parallel flow damping rate ν_∥.²⁶ The $Γi, Π∥$⁠, and $ΠE$ terms respectively capture the divergences of ion density, parallel/toroidal momentum, and $E×B$/poloidal momentum fluxes due to unresolved Fourier components. The $∂x〈ΠE〉$ form follows from Eq. (3), with $ΠE≈mini0uEiθuEix$ plus finite-Larmor-radius (FLR) corrections.¹² Linearizing Eqs. (5)–(7) and neglecting $Γi, Π∥$⁠, and $ΠE$ yields a simple dispersion relation, $ω2−2Ta/miR02=[ω/(ω+iν∥)](Ta/miq2R02)$⁠. For $q≫1$⁠, this dispersion contains a pair of weakly damped high-frequency GAMs $ω≈±(2Ta/miR02)1/2−iν∥/4q2$⁠. However, for $q$ near 1, as is typical in tokamak core plasmas, the GAMs damp at a significant fraction of $ν∥∼vti/qR$⁠, as seen in more detailed kinetic calculations.²⁷ 

To evaluate and physically understand the resulting toroidal momentum flux, we examine the free-energy balance for $nis, u∥c$⁠, and $uEz$⁠:

with energies $Eis≐Ta(nis)2/ni0, E∥c≐mini0(u∥c)2$⁠, and $EEz≐12ni0mi(uEz)2$⁠. Turbulence simulations show that the Reynolds stress (⁠$ΠE$⁠) typically acts as a source for $EEz$⁠, the geodesic transfer term (⁠$∝nisuEz$⁠) moves free energy from $EEz$ to $Eis$⁠, and both parallel flow excitation (⁠$∝nisu∥c$⁠) and the turbulent density flux sideband (⁠$∝∂x〈Γi sin θ〉$⁠) move energy out of $Eis$⁠.^25,28 Note next that the electron adiabatic response combined with nonzero $nis$ implies a radial $E×B$ drift $〈uEix cos θ〉=−cbT〈ϕ sin θ〉/Br=−cbTTe0nis/eni0Br$⁠. Recalling Eq. (4), this beats with $u∥c$ to cause a contribution $Πζ(G)=−2(ZTe0/Ti0)(ρi/r)miR0vtinisu∥c$ to the toroidal angular momentum flux $〈Πζ〉$⁠. Since $Πζ(G)$ is directly proportional to the parallel flow excitation term in Eqs. (8) and (9), we conclude that energy transfer from the pressure sideband to the parallel ion flow necessarily implies an outflux of co-current toroidal angular momentum. Indeed, since the turbulent flux term $−2u∥c∂x〈Π∥ cos θ〉$ will typically transfer energy out of $E∥c$⁠, we may use the statistical average $⋯¯$ of Eq. (9), $−bpTanisu∥c¯/r≳ν∥mini0(u∥c)2¯$⁠, to conservatively estimate the (signed) momentum flux as $Πζ(G)∼2(ν∥/bpΩci)(ZTe0/Ta)ni0miR0(u∥c)2¯$⁠. The fixed ratio of $Πζ(G)$ to the corresponding free-energy transfer term guarantees robust symmetry breaking whenever there is strong GAM damping acting via ion parallel flows, contributing to the total residual stress and rotation peaking that we will derive later.

What is happening here physically? First, Reynolds stress excites a poloidal $E×B$ flow. Poloidal variation of $B$ causes a divergence in the $E×B$ velocity, resulting in up-down-asymmetric density fluctuations, like those sketched in Fig. 1. The resulting poloidal electric field (due to adiabatic electron response) and ion pressure gradient jointly excite an ion flow along B. The net energy flow from Reynolds stress drive to damping via the ion parallel flow implies a positive correlation of the poloidal electric field and poloidal ion flow. The poloidal electric field also causes a weak radial $E×B$ drift. Due to the pitch of the magnetic field, the poloidal ion flow along the field corresponds to co- (counter-)current toroidal flow where the $E×B$ drift points radially outward (inward), which causes counter-current rotation peaking.

Energy transfer from nonaxisymmetric potential fluctuations to ion parallel flows can drive an even stronger toroidal momentum flux, but to understand its origin we must first discuss the field-aligned coordinates used in most gyrokinetic formulations: Consider an axisymmetric geometry with good nested flux surfaces, but otherwise arbitrary. Radial position is given by a flux-surface label $ρ$⁠, which is axisymmetric and satisfies $b̂·∇ρ=0$⁠. Poloidal position is specified by a distended but axisymmetric poloidal angle label $ϑ$⁠. The third coordinate $ξ$ is chosen so that $b̂·∇ξ=0$⁠, so it labels the perpendicular position within the flux surface. These choices are not arbitrary: The definition of $ξ$ implies that $b̂·∇=(b̂·∇ϑ)∂ϑ$ so $∂ϑ|ρ,ξ=(b̂·∇ϑ)−1b̂·∇$ contains only slow variation. The use of an axisymmetric $ρ$ and $ϑ$ implies that the partial $∂ξ|ρ,ϑ$ is proportional to a simple toroidal derivative $ζ̂·∇$⁠, since holding $ρ$ and $ϑ$ fixed is equivalent to holding $R$ and $z$ fixed. This property has two important implications. First, appropriate choice of $ξ$ allows toroidal periodicity to imply simple periodicity in $ξ$⁠. Second, $∂ξ∝ζ̂·∇$ vanishes for any axisymmetric quantity, in particular for equilibrium plasma parameters and the magnetic geometry. These properties allow one to construct symmetry arguments that the dominant toroidal angular momentum flux, due to the $∂ξϕ$ portion of $uEi$⁠, must vanish in the statistical average for leading-order local gyrokinetic formulations with up-down symmetric magnetic geometry.¹⁴ In contrast, the $∂ϑϕ$ portion of $uEi$ is unrestricted by the symmetry arguments^14,21 and indeed must break symmetry in the common case of energy transfer from $ϕ$ to ion parallel flows, as we now show.

We begin with the contribution of the higher-order part of the $E×B$ drift, in a simple, geometric way. Defining the radial and poloidal directions $ρ̂≐(∇ρ)/|∇ρ|$ and $p̂≐ζ̂×ρ̂$⁠, decompose $b̂=bTζ̂+bpp̂$⁠. Since $ρ̂×b̂=(ζ̂−bTb̂)/bp$⁠, the radial component of the $E×B$ drift is

The first term is the leading-order contribution, restricted by symmetry. Since the leading-order term is not purely perpendicular, the second (correction) term must be added to obtain the true $∇⊥ϕG$⁠. Although nominally smaller than the first term by $k∥/k⊥bp$⁠, it breaks the symmetry,²¹ as we will see in its contribution to $Πζ$⁠:

For emphasis, $Πζ(2)$ does not represent any effect of parallel acceleration, it is simply the advection of the parallel portion of toroidal angular momentum by a small but robustly symmetry-breaking portion of the $E×B$ drift.

To understand the momentum flux caused by $Πζ(2)$ we must examine the free-energy balance for Eqs. (1)–(3):²² 

with fluctuating pressure $Ens≐∫dV 12(Ts0/ns0)ns2$ and parallel flow $E∥s≐∫dV 12msns0u∥s2$ free energies, $E×B$ energy $EE≐∫dV ∑s12(ns0/Ts0)Z2e2ϕ(1−Γ0s)ϕ$ (includes FLR corrections), and volume integral $∫dV$⁠. Boundary terms have been assumed to vanish. The key point here is that the momentum flux term $Πζ(2)$ is directly proportional to the electrostatic acceleration of ion parallel flows, $Tϕi∥≐−Zeni0u∥i∇∥ϕG$ in Eq. (14). Although this term is sometimes referred to as ion Landau damping, it is in fact conservative, representing a transfer of energy from the potential $EE$ to ion parallel flow $E∥i$⁠. In cases with damping of turbulence via parallel ion flows, this term will tend to be positive, so that $Πζ(2)$ transports co-current momentum outward, corresponding to counter-current rotation peaking,²⁹ c. f. Fig. 2. This will especially occur when there are density fluctuations at low $k⊥$⁠, due to low frequencies, electron adiabatic response and low-$k⊥$ quasineutrality $ene0ϕ̃≈Te0ñe≈ZTe0ñi$⁠, which reduces $−∫dV u∥i[Ti0∇∥ni+Zeni0∇∥ϕG]→−Zeni0(1+Ti0/ZTe0)∫dV u∥i∇∥ϕ≈mini0∫dV u∥iD∥i≥0$⁠. Since the ion parallel flows are excited predominantly at low $k⊥$⁠, we may use this with $D∥i→ν∥u∥i$ to estimate $∫dV Πζ(2)≈(ZTe0/Ta)(ν∥/Ωciθ)mini0R0∫dV u∥i2$⁠, in which $ν∥/Ωciθ∼ρiθ/qR0$ for poloidal ion gyroradius $ρiθ≐vti/|Ωciθ|$⁠, (signed) $Ωciθ≐ZeBp/mic$⁠, and $ν∥∼vti/qR0$⁠. Alternatively, if a fraction $fL≐[∫dV Tϕi∥]/[−∑sTs0∫dV nsuEs·∇lnns0]$ of the free energy passes through $Tϕi∥$⁠, one obtains

with pressure gradient scale lengths $Lpe=Lpi≐Lp$ and turbulent radial electron $Qe$ and ion $Qi$ heat fluxes.³⁰ Defining the Prandtl number $Pr≐χφ/χi$ and assuming $χe=χi$⁠, this drives a gradient in the rigid rotation $Uζ$⁠:

Analogous manipulations in a gyrokinetic formulation lead to the same result: energy transfer from $EE$ to $E∥i$ necessarily implies an exhaust of co-current momentum, with the same basic properties, magnitude, and scaling.²³ 

A few comments: $Πζ(2)$ is a residual stress, following from symmetry-breaking due to energy transfer from $EE$ to $E∥i$⁠, regardless of the background rotation profile. The symmetry breaking is statistical: it occurs simply because free energy flows through phase space from sources to sinks. In particular, $Πζ(2)$ is nonlinear, not quasilinear—it follows from the energy transfer term summed over all modes (including damped ones such the GAMs)³¹ and does not depend on the linear mode structure of any particular instability. However, it does require the presence of fluctuations (unstable or damped) at low enough frequency to excite ion parallel flows: $ω≲vti/qR$ allows $fL∼O(1), ω≫vti/qR$ implies $fL≪1$⁠. It survives in a radially local (fluxtube) limit, not requiring any radially global effects. Although $Πζ(2)$ results from a higher-order part of the $E×B$ drift (⁠$uEi(2)$⁠), which should have little direct impact on the leading-order turbulence, it is slaved—or “parasitic”—to free-energy fluxes that are determined entirely by the leading-order physics. It can therefore be estimated even by simulations that neglect $uEi(2)$⁠, simply by evaluating $Tϕi∥$ and $fL$ a posteriori.

Although quantitative evaluation requires nonlinear simulation, we may qualitatively compare $Πζ(2)$ with experimental rotation observations. Assuming a flat current profile and letting $a$ denote device minor radius, the scaling for the mid-radius velocity gradient is roughly

which resembles Rice scaling (⁠$∝1/Ip$⁠) and has a magnitude comparable with experimental observations [strongly hollow profiles may be explained by $fL∼O(0.1)$].^5–7 Also, in ASDEX-Upgrade (AUG), counter-current momentum peaking has correlated strongly with density peaking across many discharge types.⁸ The relation may be more coincidental than causal: density peaking tends to occur due to electron precessional resonance for fluctuations with $ω≲vti/R$⁠,³² which (at core $q∼1$⁠) are the same modes that can excite ion parallel flows, thus driving counter-current peaking. Interestingly, on Alcator C-mod, the presence of counter-current peaking is correlated with the disappearance of broadband high-$k⊥$ density fluctuations.⁷ Viewed theoretically, dominant dissipation via low-$k⊥$ ion parallel flows, which implies counter-current rotation peaking in the present model, would also imply the reduction or elimination of a strong direct cascade of density fluctuations to high $k⊥$⁠, consistent with C-mod measurements. Further qualitative and quantitative comparisons are needed.

In conclusion, a geometrical correction to the $E×B$ drift causes an outward flux of co-current momentum whenever electrostatic potential energy is transferred to ion parallel flows. The robust symmetry breaking follows from the free-energy flow in phase space and does not depend on assumed linear eigenmode structure. The resulting rotation peaking is counter-current and scales with temperature over $Ip$ [Eqs. (17) and (18)]. This peaking mechanism can only act when fluctuations are low-frequency enough to excite ion parallel flows, which may explain some recent experimental observations.^7,8

Helpful discussions with C. Angioni, E. Fable, G. Hammett, P. Helander, J. Krommes, and B. Scott, and funding by the Max-Planck/Princeton Center for Plasma Physics and U.S. Dept. of Energy, Office of Science, Contract No. DE-AC02-09CH11466 are gratefully acknowledged.