The influence of sonic toroidal rotation on gyrokinetic stability and transport is studied, with important implications for heavy impurity dynamics. When centrifugal drifts and electrostatic trapping corrections are included, significant modifications to the calculated transport of heavy impurities are observed. These high-rotation corrections add to the standard Coriolis drift and toroidal rotation shear drive which are normally included in gyrokinetics. Yet, because of their complexity, centrifugal and electrostatic trapping terms (quadratic in the main ion Mach number) are not generally included in gyrokinetic codes. In this work, we explore the implications of using reduced descriptions of the rotational physics. For heavy impurities such as tungsten, cross terms due to the centrifugal force can dominate the rotation dynamics, and neglecting them is shown to lead to large errors in the impurity particle flux.

## I. INTRODUCTION

External torque driven by neutral beam injection in tokamak plasmas can drive plasma flows comparable to the ion sound speed. The presence of this toroidal plasma flow can have a profound effect on the intensity of drift-wave turbulence and corresponding levels of radial transport. In particular, the radial variation (or shear) in the **E **×** B** velocity suppresses turbulence by acting to twist plasma eddies.^{1,2} In contrast, toroidal flow shear provides a Kelvin-Helmholtz-type linear drive that can increase turbulent fluctuation levels.^{3} It is common in gyrokinetics to consider the *weak rotation* limit,^{4,5} retaining only the **E **×** B** flow, Coriolis drift, and toroidal rotation shear. In the gyrokinetic equations, these effects are represented by terms linear in the Mach number, and they provide the dominant symmetry-breaking mechanism for momentum transport. However, correct treatment of the *sonic rotation* regime requires the additional consideration of centrifugal effects.^{6–9} Sonic toroidal flow, on the order of the thermal speed, produces a centrifugal force that pushes ions toroidally outward, causing them to redistribute non-uniformly around a flux surface. As a result of quasi-neutrality, a poloidally varying electrostatic potential is generated to balance the density asymmetry, pulling the electrons with the ions.^{10} The centrifugal force induces a centrifugal drift, while the poloidally varying electrostatic potential produces modifications to the **E **×** B** drift and to the particle trapping. Because of their complexity, these new sonic terms (quadratic in the main ion Mach number) are generally ignored in gyrokinetic codes. In this work, we will show that they are nevertheless essential for the study of transport of heavy impurities, such as tungsten from the divertor wall, as the influence of the centrifugal force is amplified by their greater mass.

The first gyrokinetic code to correctly implement sonic rotation was GKW.^{11} With GKW, it has been shown that the centrifugal force leads to an enhanced mirror force, which increases the trapped electron drive and reduces the residual zonal flows, yielding a small increase in the ion heat diffusivity for ion-temperature-gradient (ITG)-dominated turbulence.^{12} However, centrifugal effects can be large for heavy impurities, such as tungsten, even when they are negligible for the bulk plasma.^{13} Studies of nonlinear turbulent transport are more limited and have focused on experimental scenarios.^{14}

In the present work, we will describe the CGYRO^{15} implementation of sonic rotation, including all quadratic Mach number terms. A major focus of the present work will be to understand the implications of neglecting centrifugal effects by comparing results from the complete sonic formulation to the simpler weak rotation approximation. We remark that the CGYRO implementation complements the rotation formulation in the drift kinetic code NEO.^{16–18} NEO implements sonic toroidal rotation using a physically consistent and mathematically equivalent ordering hierarchy.^{8,19} That is, the equations solved in CGYRO complement those solved in NEO insofar as together they represent the complete first-order deviation of the plasma from a local Maxwellian.

The rest of this paper is organized as follows. In Sec. II, the basic simulation model for full sonic rotation is described. Simulation results exploring the impact of rotation and rotation shear on gyrokinetic linear stability are shown in Sec. III. Benchmarks with GKW are also presented for verification purposes. In Sec. IV, the impact of rotation on heavy impurity transport, such as tungsten and carbon, is explored with both quasilinear and nonlinear simulations. Finally, a brief summary is given in Sec. V.

## II. HIERARCHY OF EQUATIONS

We adopt the non-orthogonal field-aligned coordinate system $(\psi ,\theta ,\alpha )$, together with the Clebsch representation for the magnetic field^{20}

The angle *α* is written in terms of the toroidal angle $\phi $ as $\alpha \u2250\phi +\nu (\psi ,\theta )$, where the function $\nu (\psi ,\theta )$ is defined with respect to the current function $I(\psi )$ by Eq. (9) in Ref. 21. In Eq. (1), *ψ* is the poloidal flux divided by $2\pi $, and *θ* simultaneously refers to an angle in the poloidal plane (at fixed $\phi $) and a parameterization of distance along a field line (at fixed *α*). In these coordinates, the Jacobian determinant and the flux-surface average are, respectively,

We consider solution of the gyrokinetic equation using the sonic rotation formalism of Sugama,^{8} which allows for flow velocities on the order of the ion thermal speed. The formalism is analogous to the gyrokinetic equations in the comoving frame of a toroidally rotating plasma developed by Peeters *et al.*^{9} on which the GKW code is based. The distribution function and fields are expanded in powers of $\rho *=\rho i/a$, the ratio of the ion gyroradius to the system size. When the rotation is sonic, the electrostatic potential must be one order larger in $\rho *$ than for the diamagnetic ordering usually assumed in kinetic theory, i.e.,

It has previously been shown by Hinton and Wong^{19} that in an axisymmetric system (considered here), $\Phi \u22121$ is a flux function, such that the zeroth-order flow velocity

is purely toroidal and species-independent. Here, $\omega 0(\psi )\u2250\u2212c\u2009d(\Phi \u22121)/d\psi $ is the angular rotation frequency.

### A. Rotating equilibrium

The $O(1)$ equilibrium equation guarantees that the zeroth-order distribution function retains the Maxwellian form, but now with respect to the rotating velocity frame

Above, $ua=v/(2vta)$ is the normalized speed in the rotating velocity frame with respect to the thermal speed $vta=Ta/ma$, such that $v2=v\u22252+v\u22a52$. While the equilibrium temperature $Ta=Ta(\psi )$ in Eq. (5) is a flux function, the equilibrium density *n _{a}* is not, varying as $na(\psi ,\theta )=Na(\psi )e\u2212\lambda a$, where

*λ*is the dimensionless effective potential energy, defined as

_{a}Here, $\Phi \u03030=\Phi 0\u2212\u27e8\Phi 0\u27e9$ is the poloidally varying part of the equilibrium potential, and we define the species Mach number as

where *R*_{0} is the midplane major radius of the flux surface. While $\u27e8\Phi 0\u27e9$ [like $Ta(\psi )$ and $Na(\psi )$] is undetermined, $\Phi \u03030$ is determined by the $O(1)$ Poisson equation

This equation is equivalent to the usual quasi-neutrality relation. In general, for an arbitrary number of species, solution of Eq. (8) for $\Phi \u03030$ requires solving a nonlinear algebraic equation. In CGYRO, we solve this numerically using Newton's method (we remark that the same equation must also be solved in NEO). For the simple case of a pure plasma, Eq. (8) can be solved analytically, yielding

The poloidally varying potential results from the density asymmetry created by the non-uniform distribution of ions around a flux surface due to the centrifugal force and pulls the electrons with the ions. Note that it is quadratic in the main ion Mach number. Regarding the input density, we find it convenient for contact with experiment to use the values at the outboard midplane. Note that, without loss of generality, we can rewrite the density at a given flux surface *ψ* as

where $\Phi *(\theta )\u2250\Phi 0(\theta )\u2212\Phi 0(\theta 0)$ and we take $\theta 0=0$. Thus, given $na(\theta 0)$, Eq. (8) can be solved for $\Phi *(\theta )$. In the diamagnetic limit (*M _{a}* = 0), $\Phi *=0$.

The magnetostatic equilibrium is also altered by the rotation. Momentum balance and Ampère's law determine the basic set of magnetostatic equilibrium equations. We can write these as

where the effective equilibrium pressure gradient is defined as

Here, $p=\u2211anaTa$ is the total equilibrium kinetic plasma pressure. Note that this is poloidally varying due to the poloidal variation of the equilibrium density. While approximations such as trace impurities and adiabatic electrons are often made regarding the total pressure, in CGYRO, the exact total pressure is used, and thus, the magnetostatic equilibrium is fully consistent with the equilibrium density and potential used in the kinetic equations. The effective pressure gradient modifies not only the gyrokinetic drift velocity (as will be detailed shortly) but also the local Grad-Shafranov equilibrium calculation.^{21} These effects may be important in the KBM-dominated plasma edge.

### B. Gyrokinetic equation

The gyrokinetic equation yields the $O(\rho *)$ fluctuating distribution in terms of the distribution function of gyrocenters, $ha(R)$. We write the equation in $(\xi ,ua)$ velocity-space coordinates which are defined with respect to the rotating velocity frame described in Sec. II A. Here, $\xi \u2250v\u2225/v$ is the cosine of the pitch angle and $ua\u2250v/(2vta)$ is the normalized velocity. For the spatial dimensions $(r,\alpha ,\theta )$, we write the distribution function $ha(R,\xi ,ua)$ in terms of the discrete indices (*p*, *n*) which are the Fourier wave numbers corresponding to the $(r,\alpha )$ coordinates, such that

where $x\u22502\pi (r\u2212r0)/L\u2208[\u20090,2\pi )$ and *L* is the radial domain size. In this notation, introducing unit vectors $b\u2250B/B$, $ex\u2250\u2207r/|\u2207r|$ and $ey\u2250b\xd7ex$, we have the corresponding wavenumbers

where $k\u22a52=kx2+ky2$. The related poloidal wavenumber $k\theta $ is defined so that it is a proper flux-surface function: $k\theta \u2250nq/r$, where $q(\psi )$ is the safety factor. The field potential $\Psi a(R)$^{15} can also be expanded as

where

Here, $\gamma a\u2250k\u22a5v\u22a5/\Omega ca$. Note that we have defined a new effective electrostatic potential as

which is just the potential in the rotating frame, using the Lorentz transformation for the electric field: $\delta Eeff=E+U0/c\xd7B$. The nonadiabatic distribution $H\u0303a$ is defined in terms of $\Psi \u0303a$ as

We write the gyrokinetic equation in dimensionless form by normalizing to the deuteron sound speed $cs=Te/mD$ (where *T _{e}* is the electron temperature and

*m*is the deuteron mass) and the midplane minor radius of the last closed flux surface,

_{D}*a*. Thus, we define the normalized time as $\tau \u2250(cs/a)t$. The dimensionless gyrokinetic equation for $h\u0303a(n,p,\theta ,\xi ,ua)$ becomes

In Eq. (21), $CabL$ is the linearized gyrophase-averaged collision operator.^{22} The spectral representation of the nonlinear term $h\u0303a\u2217\Psi \u0303a$ (a generalized convolution) is given by Eq. (56) of Ref. 15. The dimensionless operators are

**streaming**

**trapping**

**drift**

**gradient drive**

In Eqs. (22)–(25), the species gyroradius is defined as $\rho a\u2250vtamac/(zaeB)$, with the effective ion-sound gyroradius given by $\rho s\u2250csmDc/(eBunit)$, where $Bunit(r)=(q/r)\psi \u2032$ is the effective magnetic field. We have also introduced an effective velocity $ws\u2250cs(J\psi B)/a$ and the density and temperature gradient inverse scale lengths, $1/Lna=\u2212d\u2009ln\u2009na(\theta 0)/dr$ and $1/LTa=\u2212dlnTa/dr$.

Note that the rotation terms have been written explicitly in Eqs. (22)–(25). We first define the *weak rotation* terms, which are linear in *M _{a}*,

**Coriolis drift**

**rotation shear drive**

Here, $\gamma p=\u2212R0d\omega 0/dr$ is the toroidal rotation shearing rate.^{23} Note that *γ _{p}* is sometimes referred to by the misnomer “parallel velocity shearing rate,” which reflects the historical usage in sheared slab geometry for which $\gamma p=\u2212R0d(U\u2225/R)/dr$.

^{24}The sheared

**E**×

**B**flow Ω

_{s}in Eq. (21) is also a weak rotation term, being linear in

*M*. For inclusion of sheared

_{a}**E**×

**B**flow in Eq. (21), we first define the dimensionless

**E**×

**B**shearing rate

where $\gamma E=\u2212(r/q)d\omega 0/dr$ is the Waltz shearing rate^{1} and *γ _{s}* is the domain shearing rate. In Eq. (21), the total shearing term is $\u2212i\Omega sXh\u0303a$, where

*X*is a spectral shearing operator, defined by Eq. (21) of Ref. 25.

For the case of *sonic rotation*, the centrifugal terms must also be considered. These terms are quadratic in *M _{a}* and contribute to the

**trapping**operator through

the **drift** operator through

and the **gradient drive** operator through

Note that centrifugal effects also modify the equilibrium density in the equilibrium distribution $f0a$ in $ha(x,\alpha )$ in Eq. (14) and in the effective pressure gradient in the drift in Eq. (24). It is common in gyrokinetics to neglect the centrifugal effects due to their complexity and to instead only include the weak rotation terms, namely, the Coriolis drift, parallel velocity shear drift, and sheared **E **×** B** flow. The implications of this reduction will be explored in detail shortly.

The dimensionless, spectral gyrocenter-Maxwell equations are given by

Here, we have introduced the Debye length $\lambda D=Te/(4\pi nee2)$ and the effective electron beta $\beta e,unit=8\pi neTe/Bunit2$. Note that with sonic rotation, we can solve the Maxwell equations for the effective electrostatic potential $\delta \varphi \u0303eff$ and thus the field potential given by Eq. (18). In the limit of weak rotation, $\delta \varphi \u0303eff\u2192\delta \varphi \u0303$.

Finally, the turbulent radial particle flux Γ_{a}, energy flux *Q _{a}*, and momentum flux Π

_{a}are defined in accordance with Ref. 8 as

where

The field quantity $\chi \u0303a$ is given by

Note that the $(n=0,p=0)$ component of the distribution, $h\u0303a(0,0,\theta )$, is zero in gyrokinetic theory. Thus, spatial averages of the distribution, such as the gyrokinetic bootstrap current, are zero.

To recap, CGYRO solves a 5 D nonlinear, electromagnetic gyrokinetic equation, Eq. (21), for each species *a*. These equations couple through the gyrocenter Maxwell equations, Eqs. (32)–(34), and also through the collision operator. As described in Ref. 15, the solver is both pseudospectral in velocity space dimensions $(\xi ,ua)$ and spectral in the radial and binormal spatial dimensions $(x,\alpha )$ and uses a 6th-order finite difference scheme in *θ* (the fieldline direction) with a 6th-order conservative dissipation algorithm.

## III. SIMULATION PARAMETERIZATION RESULTS

The simulations use the General Atomics (GA) standard case parameters,^{26} considering a single ion species (deuterium) and electrons, with $mi/me=60$ with the following: $R0/a=3,\u2009r/a=0.5$, *q* = 2, *s* = 1, *T _{i}* =

*T*, $a/Lni=a/Lne=1,\u2009a/LTi=a/LTe=3$. In CGYRO, we specify the density and the density gradient at the outboard midplane, $\theta 0=0$. For the geometry, we use an unshifted Miller local circular equilibrium.

_{e}^{27}We emphasize that the local equilibrium is not a model equilibrium, but rather an exact solution of the Grad-Shafranov equations on the specified flux-surface. The local equilibrium method is complicated to implement and the reader should consult Ref. 21.

To simulate a single linear eigenmode, we consider fixed $k\theta \rho s=0.3$. To vary the rotation frequency *ω*_{0} and rotation shear $d\omega 0/dr$, we (equivalently) specify the main ion Mach number $Mi=\omega 0R0/vti$ and the toroidal rotation shearing rate $\gamma p=\u2212R0d\omega 0/dr$. The species-dependent Mach numbers *M _{a}* are varied consistently with

*M*. The poloidally varying equilibrium density, computed by CGYRO, is shown in Fig. 1 for varying

_{i}*M*. Accumulation of ions at the outboard midplane (

_{i}*θ*= 0) is observed as

*M*increases due to the centrifugal force. To assess the impact of rotation, three cases are considered:

_{i}**Sonic rotation:**All sonic rotation terms are included in the equilibrium and gyrokinetic equations.**Weak rotation:**Only rotation terms which are linear in*M*are included in the gyrokinetic equations; specifically, the Coriolis drift $\Omega d,cor$, the toroidal rotation shear drive $\Omega *,trs$, and the_{a}**E**×**B**term Ω_{s}. No centrifugal terms are included, meaning $\Phi *=0$ and $na=na(0)$.**Only centrifugal terms:**Only rotation terms due to the centrifugal force are included in the gyrokinetic equations; specifically, $\Omega \xi ,cf,\u2009\Omega d,cf$, and $\Omega *,cf$. The poloidal variation of the equilibrium density is included in $f0a$ as well as in the pressure-gradient drift $(\u2207p)eff$. Weak rotation terms are neglected.

As we have explained previously, the limit of weak rotation is commonly assumed in linear and nonlinear gyrokinetic simulations. It is equivalent to retaining only terms linear in the Mach number *M _{a}*. While this approximation is typically very accurate for deuterium ions at low rotation, it is not generally valid for heavy impurities.

### A. Benchmarks with GKW

To test the implementation of all sonic terms, CGYRO has been benchmarked with GKW^{11} using the GA standard case parameters specified earlier. The results are summarized in Figs. 2 and 3. The results show that the intercode agreement is excellent, without any significant deviation, in scans over *M _{i}* and

*γ*in both low- and high-rotation regimes.

_{p}### B. ITG and trapped electron mode (TEM) linear stability

The effect of rotation on gyrokinetic linear stability of ion-temperature-gradient (ITG) modes and trapped electron modes (TEM) using the GA standard case parameters is further studied in Fig. 4. The results with sonic rotation are compared with those including only weak rotation terms and those including only centrifugal terms. At low *M _{i}*, the ITG mode is mostly affected by the stabilizing Coriolis drift $\Omega d,cor$, while centrifugal terms are insignificant (i.e., red curve follows the blue curve). However, as

*M*increases, the ITG mode transitions to a TEM, which is further destabilized as

_{i}*M*increases. For the TEM, the centrifugal terms (green curve) dominate the Coriolis drift term (blue curve). Although the centrifugal force itself is small for the electrons due to their small mass, enhancement of the mirror force due to the induced poloidally varying electrostatic potential leads to an increased effective electron trapped fraction. For these scans, a resolution of ($N\theta =N\xi =24$) is required to resolve the enhanced particle trapping. This trapping effect has a strong destabilizing influence on the TEM, causing the mode transition to occur. This has been previously shown with the GKW code for similar parameters.

_{i}^{12}However, the effect on TEMs can be quickly reduced by the detrapping of electrons due to collisions, as shown in Fig. 5. The results use the Sugama collision model. As the collision frequency increases, the growth rate of the TEM decreases and eventually mode transition to a dominant ITG mode occurs.

In contrast to the effect of the main ion Mach number, Fig. 6 shows that increasing the toroidal rotation shear *γ _{P}* destabilizes the ITG mode. Comparison with the

*M*= 0 case indicates that, at large

_{i}*γ*, the $Mi\gamma p$ cross term explicitly in $\Omega *,cf$ and implicitly in $\u2207\Phi *$ in $\Omega d,cf$ enhances the destabilizing effect. In contrast to the ITG mode, the TEM is stabilized by the toroidal rotation shear, leading to TEM to ITG mode transition.

_{p}## IV. TURBULENT TRANSPORT OF HEAVY IMPURITIES

To study the effects of rotation on turbulent transport, we consider the GA standard case parameters used above with the addition of a trace impurity species. For reactors, transport of heavy impurities from the material surface into the high temperature core must be minimized since accumulation can lead to fuel dilution and radiation losses. Both fully ionized carbon (*z _{i}* = 6) and partially ionized tungsten (

*z*= 41) are considered. The trace limit is assumed, with $nimp/ne=10\u22125$. The impurity species is assumed to be in thermal equilibrium with the main ions. We consider a typical fixed rotation frequency of $Mi=0.2$ and varying rotation shear. The results compare inclusion of full sonic rotation effects with those including only weak rotation terms, neglecting centrifugal effects. The sonic and weak rotation cases converge when

_{i}*M*= 0. The poloidal variation of the equilibrium density for the ions in the full sonic rotation case is shown in Fig. 7. Despite the small rotation frequency, the density asymmetry for the impurities is large, especially for tungsten. This is not surprising since $MC=0.49$ and $MW=1.92$.

_{a}### A. Quasilinear critical gradient

The results for the critical impurity density gradient scale length are shown in Fig. 8, comparing carbon and tungsten. The critical impurity density gradient length scale is defined as that for which the quasilinear turbulent radial particle flux is zero. For carbon, the case with only weak rotation terms (blue curve) qualitatively follows the trend of the full rotation case (red curve), with the latter predicting a larger critical gradient. The dominant rotation effect is the toroidal rotation shear drive $\Omega *,trs$. The critical impurity density gradient decreases as *γ _{p}* increases, such that the impurity density profile is peaked for small rotation shear, $\gamma pa/cs<4$, and hollow for large rotation shear. In contrast to carbon, for tungsten, including only weak rotation terms while neglecting the centrifugal terms produces incorrect results. This is not surprising considering the large poloidal density asymmetry due to the centrifugal force shown in Fig. 7. For $\gamma pa/cs>2$, the case with inclusion of only weak rotation terms (blue curve) predicts a hollow impurity density profile, similar to the case with carbon for large

*γ*. However, inclusion of sonic rotation (red curve) yields a critical impurity density profile which is peaked for the entire range of

_{p}*γ*. The fact that both the weak and sonic rotation curves are close to the

_{p}*M*= 0 limit (green curve) at $\gamma p=0$ indicates that the main effect is due to the $Ma\gamma P$ cross-terms rather than the $Ma2$ terms.

_{a}### B. Nonlinear turbulent flux

Figures 9 and 10 show the corresponding nonlinear turbulent fluxes for $a/LnW=0$. The particle and energy fluxes are normalized by the gyroBohm factors

and by the relative density, such that ${\Gamma a,Qa}norm={\Gamma a,Qa}/(na/ne)$. The deuterium ions and electrons are negligibly affected by the centrifugal force, and thus, the full sonic and weak rotation results are very similar for these species. However, similar to the quasilinear results, the tungsten particle flux is strongly affected by the centrifugal dynamics and thus including only weak rotation terms gives a large error. With full rotation, the tungsten particle flux is radially inward, with its magnitude increasing strongly with increasing rotation shear. This indicates a strong inward pinch, which leads to detrimental core accumulation in a reactor. Including only weak rotation terms incorrectly predicts an outward particle flux at large rotation shear, as shown by the blue curve in Fig. 10(b). Note that inclusion of **E **×** B** shear (not shown here) has a negligible effect on these results. Figure 10(a) shows that the tungsten energy flux is destabilized by the rotation shear, in contrast with the slight stabilizing effect observed for the deuterium ions and electrons in Fig. 9. The insignificance of centrifugal dynamics indicates that, unlike the tungsten particle flux, the tungsten energy flux is dominated by the toroidal rotation shear drive $\Omega *,trs$. For comparison, the results for the carbon nonlinear turbulent fluxes for $a/LnC=0$ are shown in Fig. 11. As expected, the carbon fluxes are less affected by the centrifugal dynamics than tungsten but more strongly affected than the deuterium ions. Similar to tungsten, an inward pinch is predicted with full sonic rotation, and neglecting centrifugal dynamics incorrectly predicts an outward flux, as shown by the blue curve in Fig. 11(b). Ultimately, these results indicate that inclusion of full sonic rotation, including centrifugal effects, is critical for heavy impurities, such as tungsten, as well as for moderate impurities like carbon at large rotation shear.

Finally, Fig. 12 compares the impurity turbulent and neoclassical fluxes. The neoclassical fluxes are computed with NEO,^{18} which also includes full sonic rotation. For low collision frequency typical in the tokamak core ($\nu \xafe=0.1$), the neoclassical fluxes are much smaller than the turbulent fluxes, by 1–2 orders of magnitude. However, since the magnitude of the neoclassical fluxes increases with increasing $\nu \xafe$ while the magnitude of the turbulent fluxes decreases, at $\nu \xafe$ = 1, which is typical in the H-mode pedestal, the neoclassical fluxes become competitive with the turbulent fluxes. In fact, for $\gamma pa/cs<2$, the neoclassical fluxes dominate. Thus, both must ultimately be considered for accurate transport modeling of tungsten. The neoclassical tungsten particle fluxes are also directed inward, enhancing the pinch effect.

## V. SUMMARY

The impact of sonic toroidal rotation on gyrokinetic stability is studied using the CGYRO code, which solves the gyrokinetic equation using the sonic rotation formalism of Sugama.^{8} This formalism is based on a rigorous asymptotic expansion of the Fokker-Planck equation and is valid for flow velocities on the order of the ion thermal speed. Centrifugal drifts and electrostatic trapping dynamics are included, in addition to the standard *weak rotation* terms included in most gyrokinetic codes: the **E **×** B** flow, Coriolis drift, and toroidal rotation shear. Because of their complexity, these sonic rotation effects—which are quadratic in the Mach number—are rarely considered in gyrokinetic simulations.

For a pure plasma, the ITG mode is primarily affected by the stabilizing Coriolis drift for low ion Mach number ($Mi<1$). At large *M _{i}*, the dominant effect of sonic rotation is on the electron dynamics. Though the centrifugal force is small for electrons due to their small mass, enhancement of the mirror force leads to an increased effective electron trapped fraction. This has a strong destabilizing influence on TEMs and leads to mode transition from ITG to TEM for increasing

*M*. However, the effect on TEMs can be quickly reduced by the detrapping of electrons due to collisions. In contrast, the toroidal rotation shear

_{i}*γ*is destabilizing for the ITG mode and stabilizing for the TEM, leading to a TEM-to-ITG mode transition. For heavy impurity turbulent transport at moderate

_{p}*M*, including only the weak rotation terms (i.e.,

_{i}**E**×

**B**flow, Coriolis drift, and toroidal rotation shear) while neglecting the centrifugal terms results in a large error. For analysis of the quasilinear critical tungsten density gradient, inclusion of only weak rotation terms predicts a hollow impurity density profile, while the complete sonic theory yields a critical impurity density profile which is peaked due to the $Ma\gamma P$ cross-terms. In contrast, carbon is not significantly affected by rotation, and a hollow impurity density profile is predicted at large rotation shear.

Nonlinear turbulent fluxes for $a/LnW=0$ similarly show that the sonic formulation must be used, yielding a strong inward pinch, which can lead to detrimental impurity accumulation in the reactor core. Including only weak rotation terms incorrectly predicts an outward particle flux at large rotation shear. At a large collision rate, the neoclassical particle fluxes become competitive with the turbulent fluxes, enhancing the pinch effect. NEO simulations in conjunction with *quasilinear* GKW calculations have previously been used to analyze hybrid scenarios in JET. In this work, it was found that the neoclassical transport rather than the turbulent transport dominates the central part of the plasma due to enhancement of the neoclassical diffusive and convective components by rotation.^{28} More recent results, however, indicate that quasilinear simulations may be inaccurate for mixed turbulence regimes with comparable ion and electron transport.^{14} Unfortunately, impurity transport modeling that couples neoclassical with nonlinear gyrokinetic simulation and incorporates heavy impurity sources and radiative losses is still computationally impractical. Thus, it is important that reduced models of transport^{29,30} be extended to treat the sonic regime correctly in order to provide realistic assessments of heavy impurity accumulation.

## ACKNOWLEDGMENTS

This material is based on work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, Theory program, under Award No. DE-FG02-95ER54309, and by the Edge Simulation Laboratory project under Grant No. DE-FC02-06ER54873. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of the authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.