Physical considerations

Sonic rotation can arise in tokamaks from torque due to neutral beam injection. It is believed to play an important role in the suppression of turbulence and the formation of transport barriers through \(\mathbf{E}\times\mathbf{B}\) shear. In addition, it produces a strong centrifugal force that pushes ions toroidally outward, causing them to redistribute nonuniformly around a flux surface. As a result of quasi-neutrality, a poloidally-varying electrostatic potential is generated by the electrons to balance the density asymmetry. This generates centrifugal drifts in addition to the usual \(\mathbf{E}\times\mathbf{B}\) rotation, Coriolis drift, and parallel velocity shear. Because of their complexity, these centrifugal effects, which are second-order in the main ion Mach number, are ignored in most neoclassical and gyrokinetic codes. Inclusion of full sonic rotation effects is particularly critical for studying transport of heavy impurities [BC08], as the influence of the centrifugal force is amplified by their large mass.

Sonic rotation formalism

The rigorous derivation of the gyrokinetic equations for sonic rotation was carried out by Sugama [SH98]. Sugama followed the underlying rotation formulation of Hinton and Wong [HW85], who showed that in an axisymmetric system the zeroth-order flow velocity is

\[\mathbf{U}_0 = \omega_0(\psi) R^2 \nabla \varphi \; .\]

It is this scalar flux-function, \(\omega_0\), that is the true free (input) function in both neoclassical and gyrokinetic theory. Here, \(\psi\) is the poloidal flux.

The rotation profile

The relevant profile that gives a complete specification of rotation is the angular frequency \(\omega_0(r)\), which in defined consistently throughout GACODE as

\[\omega_0(r) \longrightarrow \frac{c E_r }{R B_p} \; .\]

Thus, if the \(E_r\) is known from a priori analysis (say, through velocity measurements and force balance), then \(\omega_0\) is given by the formula above. This equality is subtle (hence the reason for the arrow rather than an equality) and is clarified in detail in the last section.

Input parameters

CGYRO and NEO implement full sonic rotation (GYRO implements only a reduced model) according to the formulation of Hinton and Wong [HW85]. Neoclassically, the induced poloidally-varying electrostatic potential leads to the formation of potential wells. In the banana regime these increase the effective trapped particle fraction, and in the Pfirsh-Schlüter regime increase the effective toroidal curvature. In both instances, this may lead to enhanced neoclassical transport. The code inputs are given in the tables below.


Note that all the inputs are derived from the single free function, \(\omega_0\).

input.cgyro parameter




\(\displaystyle \frac{a}{c_s} \, \gamma_{\rm E} \doteq -\frac{r}{q}\frac{d \omega_{0}}{d r}\)

\(\mathbf{E}\times\mathbf{B}\) shearing rate


\(\displaystyle \frac{a}{c_s} \, \gamma_p \doteq -R_0\frac{d \omega_{0}}{d r}\)

rotation shearing rate


\(\displaystyle M \doteq \frac{\omega_0 R_0}{c_s}\)

rotation rate

input.neo parameter




\(\displaystyle \frac{a}{v_{norm}} \omega_0\)

rotation rate


\(\displaystyle \frac{a^{2}}{v_{norm}} \frac{d \omega_{0}}{dr}\)

derivative of rotation rate

Theoretical basis for sonic rotation

In presence of rapid rotation, where the flow speed \(U\) is allowed to be of the order of the ion thermal speed, the Lorentz force term in the the Fokker-Planck equation (see [HW85])

\[\frac{\partial f_i}{\partial t} + \mathbf{v} \cdot \nabla f_i + \frac{e}{m_i}(\mathbf{E}+\mathbf{v} \times \mathbf{B}) \cdot \frac{\partial f_i}{\partial \mathbf{v}} = C_i + S_i\]

becomes the leading term. Under these circumstances Hinton and Wong show that

\[\mathbf{E}_{-1} + \frac{\mathbf{U}_0}{c} \times \mathbf{B} = 0 \; ,\]

where the subscripts represents the order with respect to the drift ordering in \(\rho_i/a\). Here \(\mathbf{U}_0\) is a purely toroidal velocity and is species independent. This ordering is applied to all fields and moments

\[\begin{split}\begin{matrix} f_i & = & & & f_{i,0} & + & f_{i,1} & + & \ldots \\ \Phi & = & \Phi_{-1} & + & \Phi_0 & + & \Phi_1 & + & \ldots \\ \mathbf{U} & = & & & \mathbf{U}_0 & + & \mathbf{U}_1 & + & \ldots \\ U_\varphi & = & & & U_{\varphi,0} & + & U_{\varphi,1} & + & \ldots \\ U_\theta & = & & & & & U_{\theta,1} & + & \ldots \end{matrix}\end{split}\]

where \(U_\varphi \doteq \mathbf{e}_\varphi \cdot \mathbf{U}\) is the toroidal velocity and \(U_\theta \doteq \mathbf{e}_\theta \cdot \mathbf{U}\) is the poloidal velocity. The leading-order sonic flow is toroidal and independent of species

\[\mathbf{U}_0 = \omega_0(\psi) R \mathbf{e}_{\varphi} \quad \text{where} \quad \omega_{0}(\psi) \doteq -c \frac{d \Phi_{-1}}{d \psi} \; .\]

It is important to note that \(\Phi_{-1}\) is a flux function, whereas higher orders are not constant on a flux surface.

Connection to experimental data

We remark that \(\omega_{0}\) is a theoretical quantity that cannot be measured in the experiment. This is similar to the observation that the experimentally-measured temperature \(T_i\) is really the sum of an equilibrium temperature and a small fluctuating temperature driven by turbulence: \(T_i = T_{i0} + T_{i1}\). If the drift ordering is valid, then we are justified in approximating the equilibrium temperature \(T_{i0}\) by the measured temperature \(T_i\). For the rotation frequency, similar considerations hold. We note that the theory shows that the potential always appears in the combination

\[\Phi_{-1} + \left\langle \Phi_0 \right\rangle \; ,\]

where an angle bracket denotes a flux-surface average. By analogy with the temperature, the rotation frequency can related to the experimentally-deduced radial electric field \(E_r\), where \(E_r = -|\nabla r| d\Phi/dr\), according to

\[\omega_0 + \omega_1 \simeq \frac{c E_r}{R B_p} \quad\text{where}\quad \omega_1 \doteq -c \frac{d \left\langle \Phi_0 \right\rangle }{d\psi} \; .\]

In practice, we can set \(\omega_1 = 0\) without loss of generality and all the rotation is contained in \(\omega_0\). Alternatively, in the diamagnetic rotation limit, we set \(\omega_0 = 0\) with the rotation contained in \(\omega_1\). The present theory works consistently in both cases. Finally, the toroidal velocities \(U_{\varphi,0} + U_{\varphi,1}\) are treated in the same way.

Consistency with force balance

In experimental analyses the radial force balance relation is often used

\[E_r = \frac{R B_p}{n_a z_a e} \frac{d p_a}{d\psi} + \frac{U_\varphi}{c} B_p - \frac{U_\theta}{c} B_t \; .\]


We emphasize that this relation is valid at long wavelength (equilibrium scales) only, and is subject to the same ordering requirements as standard neoclassical and gyrokinetic theory. This means a restriction on the steepness of gradients in the form \(d \ln p/dr \ll 1/\rho_i\). See, for example, the discussion in [SWNN11].

The force balance relation contains terms of order 0 and 1, as described in the previous sections. We can write the velocities in terms of the neoclassical flow coefficient \(K_a\) (see [BC09]) as

\[\begin{split}\begin{align} U_\varphi = &~ \frac{K_a}{n_a} B_t + \omega_{1,a} R + \omega_0 R \; , \\ U_\theta = &~ \frac{K_a}{n_a} B_p \; . \end{align}\end{split}\]

In the expression for \(U_\varphi\), we have defined the angular frequencies

\[\begin{split}\begin{align} \omega_{1,a} = &~ -c \frac{d \left\langle\Phi_0\right\rangle}{d\psi} - \frac{c}{n_a z_a e} \frac{d p_a}{d\psi} + {\cal O}(M^2) \; , \\ \omega_0 = &~ -c \frac{d \Phi_{-1} }{d\psi} \; . \end{align}\end{split}\]

Substitution of the neoclassical flows into the force balance relation shows that all species-dependent terms cancel, leaving

\[E_r = \frac{R B_p}{c} \left( \omega_0 + \omega_1 \right) \; ,\]

where the species-independent frequency \(\omega_1\) is discussed in the previous section.