1 Early history

The term gyrokinetics was first used by Catto & Tsang (Reference Catto and Tsang1977) slightly before the gyrokinetic change of variables was introduced by Catto (Reference Catto1978). Catto & Tsang (Reference Catto and Tsang1977) derived a linearized collision operator for perpendicular wavelengths comparable to a gyroradius by using a radial WKB (Wentzel–Kramer–Brillouin) or eikonal approximation for the waves. Their procedure made use of the electrostatic efforts of Rutherford & Frieman (Reference Rutherford and Frieman1968) and Taylor & Hastie (Reference Taylor and Hastie1968) that treated small perpendicular wavelengths by using an eikonal technique rather than employing a gyrokinetic change of variables. However, the Catto (Reference Catto1978) gyrokinetic change of variables streamlined these earlier treatments of magnetized plasmas and was adopted and appreciated even in the era before codes (BC) since it avoided the need for complicated trajectory integral solutions of the kinetic equation. Indeed, the new gyrokinetic variables retained gyration and drift effects while removing fast time variation associated with gyrofrequency effects. Moreover, the gyrokinetic equation reduced to the Hazeltine (Reference Hazeltine1973) drift kinetic equation in the long wavelength limit.

The gyrokinetic change of variables removes the gyrophase as a velocity space variable by averaging out the rapid gyration motion associated with gyrofrequency time scales. It allows elongated structure along magnetic field lines for low frequency fluctuations as observed in magnetically confined and space plasmas, while allowing short, gyroradius wavelengths across the magnetic field lines. The remaining two velocity space variables enter as the magnetic moment adiabatic invariant and a near constant of the motion energy for typical analytic calculations, while turbulence simulations normally employ a parallel velocity variable along with the magnetic moment. The gyrokinetic orderings are consistent with equilibrium descriptions, low frequency electromagnetic stability, and the evaluation of collisional fluxes, as well as turbulent transport. Recent advances are allowing some simulations to be performed on the much longer turbulent and collisional transport time scales needed to evolve radial profiles with momentum transport retained to properly determine the axisymmetric radial electric field in a tokamak.

An early electromagnetic gyrokinetic treatment neglecting compressional effects appeared shortly after gyrokinetic variables were introduced (Hitchcock & Hazeltine Reference Hitchcock and Hazeltine1978), and soon after fully electromagnetic linear descriptions appeared, both without the use of gyrokinetic variables (Antonsen & Lane Reference Antonsen and Lane1980) and with the use of gyrokinetic variables (Catto, Tang & Baldwin Reference Catto, Tang and Baldwin1981). In addition, nonlinear electromagnetic treatments appeared that employed unperturbed gyrokinetic variables in tokamak geometry (Hitchcock & Hazeltine Reference Hitchcock and Hazeltine1978; Frieman & Chen Reference Frieman and Chen1982), and nonlinear electrostatic formulations in slab geometry (Lee Reference Lee1983) retaining perturbed gyrokinetic variables (Dubin et al. Reference Dubin, Krommes, Oberman and Lee1983). Hahm, Lee & Brizard (Reference Hahm, Lee and Brizard1988) extended the perturbed gyrokinetic variables to nonlinear electromagnetics in a uniform magnetic field $\boldsymbol{B}$ , and Hahm (Reference Hahm1988) and Brizard (Reference Brizard1989) made nonlinear generalizations to tokamak geometry.

Two early codes appeared in the 1980s as well. The first was the linear continuum eigenvalue gyrokinetic code FULL (Rewoldt, Tang & Chance Reference Rewoldt, Tang and Chance1982), and the second was the nonlinear gyrokinetic particle-in-cell (PIC) code of Lee (Reference Lee1983) in slab geometry. Years later, Kotschenreuther, Rewoldt & Tang (Reference Kotschenreuther, Rewoldt and Tang1995) built the first initial value, fully implicit linear delta f gyrokinetic code in tokamak geometry GKS (for GyroKinetic Stability), with delta f indicating that only the correction to a Maxwellian was being found gyrokinetically. GKS was extended by Dorland et al. (Reference Dorland, Jenko, Kotschenreuther and Rogers2000) to become the first electromagnetic nonlinear continuum delta f gyrokinetic tokamak turbulence code GS2. GS2 also retained trapped and passing particles and was fully parallelized.

The earliest nonlinear delta f flux tube PIC codes were electrostatic (Sydora Reference Sydora1995; Dimits et al. Reference Dimits, Williams, Byers and Cohen1996). The later work developed the code that became PG3EQ and discovered an important discrepancy between gyrokinetic and gyrofluid descriptions of ion temperature gradient (ITG) turbulence. It revealed the key role of zonal flow, which is a robust, temporally varying $\boldsymbol{E}\times \boldsymbol{B}$ or electric field $\boldsymbol{E}$ drift due to an axisymmetric radial electric field with strong radial variation, but no poloidal variation. It became apparent that zonal flow acts to regulate and limit ITG turbulent transport, especially near marginal stability. This insight led to the understanding that short wavelength turbulence gives rise to nonlinear beating that deposits some energy into the axisymmetric radial electric field variation whose poloidal variation damps to a residual steady state level to become the zonal flow. Further insights into zonal flow behaviour (absent from gyrofluid codes of the time) came from code simulations and comparisons (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000), leading to an understanding that there is a nonlinear Dimits shift away from the ITG linear stability threshold. The Dimits shift occurs because the generated zonal flow nonlinearly suppresses ITG turbulence near marginality. As a result, the onset of large turbulent transport is delayed and shifted to a nonlinear threshold value higher than the linear prediction.

An analytic linear initial value check was developed by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) to verify whether the ion polarization effects associated with magnetic drift are being properly retained in gyrokinetic codes. It yields a non-zero residual level of radial electric field (and thereby zonal flow) when the electric field is linearly perturbed. Later, Sugama & Watanabe (Reference Sugama and Watanabe2006) successfully described the transient damping of poloidal variation by geodesic acoustic modes (GAMs) to a residual steady state zonal flow. Other local delta f gyrokinetic continuum codes soon became available, like GYRO (Candy & Waltz Reference Candy and Waltz2003), GENE (Dannert & Jenko Reference Dannert and Jenko2005), GKV (Sugama & Watanabe Reference Sugama and Watanabe2006) and GKW (Peeters et al. Reference Peeters, Strintzi, Camenen, Angioni, Casson, Hornsby and Snodin2009), and even global PIC codes such as ELMFIRE (Heikkinen et al. Reference Heikkinen, Kiviniemi, Kurki-Suonio, Peeters and Sipilä2001). Additional global delta f PIC codes like GTC (Lin et al. Reference Lin, Ethier, Hahm and Tang2002), GEM (Chen & Parker Reference Chen and Parker2003), GTS (Wang et al. Reference Wang, Lin, Tang, Lee, Ethier, Lewandowski, Rewoldt, Hahm and Manickam2006) and ORB5 (Jolliet et al. Reference Jolliet, Bottino, Angelino, Hatzky, Tran, McMillan, Sauter, Appert, Idomura and Villard2007), were developed soon after GS2, PG3EQ and ELMFIRE.

2 Recent theoretical insights

Intrinsic ambipolarity means that the radial electric field cannot be determined from $\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}\rangle =0$ , or, upon integrating, from $\langle \boldsymbol{J}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle =0$ , since this constraint is automatically or intrinsically satisfied within a flux function in the electrostatic potential that is determined by momentum transport. Here $\boldsymbol{J}$ is the current density and $\langle \cdots \rangle$ denotes a flux surface average with $\unicode[STIX]{x1D713}$ the poloidal flux function. Sugama et al. (Reference Sugama, Okamoto, Horton and Wakatani1996) demonstrated that neoclassical and turbulent transport are separately intrinsically ambipolar in a tokamak. Although their treatment is fully electromagnetic, some details relied on using a quasilinear treatment and Krook collision operator. These limitations where removed electrostatically by Parra & Catto (Reference Parra and Catto2009). Appendix A of Sugama & Horton (Reference Sugama and Horton1998) also showed intrinsic ambipolarity electromagnetically in a tokamak with sonic toroidal flow within a quasilinear treatment of the Onsager symmetry. Intrinsic ambipolarity implies that the radial electric field is set by the need to satisfy toroidal angular momentum conservation in an turbulent tokamak as suggested by Catto et al. (Reference Catto, Simakov, Parra and Kagan2008) and verified by Parra & Catto (Reference Parra and Catto2008), as well as the need to satisfy important symmetry properties of the gyrokinetic equation (Parra, Barnes & Catto Reference Parra, Barnes and Catto2011; Sugama et al. Reference Sugama, Watanabe, Nunami and Nishimura2011).

Until 2008 it was widely assumed that the gyrokinetic equation and quasineutrality were all that was required to obtain a complete electrostatic solution. However, by much more carefully verifying the suggestion of Catto et al. (Reference Catto, Simakov, Parra and Kagan2008), Parra & Catto (Reference Parra and Catto2008) demonstrated that much greater care was required to describe momentum transport in a tokamak, and that very small errors would lead to fake torques (Parra & Catto Reference Parra and Catto2010b ,Reference Parra and Catto c ) and result in unphysical momentum transport. The controversy that ensued seems to have subsided as the result of more general gyrokinetic descriptions of turbulent momentum transport (Parra & Catto Reference Parra and Catto2010a ,Reference Parra and Catto b ; Parra et al. Reference Parra, Barnes and Catto2011; Parra et al. Reference Parra, Barnes, Calvo and Catto2012; Parra & Barnes Reference Parra and Barnes2015a ,Reference Parra and Barnes b ). However, the numerical implementation of a completely general description of momentum transport to the requisite order remains incomplete. The remaining subtlety is that delta f descriptions rely on scale separation, and therefore, require certain global features in the vicinity of a flux surface (Parra & Barnes Reference Parra and Barnes2015a ). Full f descriptions avoid scale separation, but cannot yet be formulated to high enough order in the gyrokinetic variables (Parra et al. Reference Parra, Barnes, Calvo and Catto2012). For example, they are roughly equivalent to trying to solve the short mean free path, long wavelength limit of the gyrokinetic equation to very high order to avoid using a Chapman–Enskog approach.

3 Kinetic equation in a tokamak

To keep the presentation as simple as possible while retaining geometric effects, the electrostatic limit of an axisymmetric tokamak is considered first. The magnetic field is taken as

(3.1) $$\begin{eqnarray}\displaystyle \boldsymbol{B}=B\boldsymbol{b}=I\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D735}(\unicode[STIX]{x1D701}-q\unicode[STIX]{x1D717})\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D713}$ , $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D717}$ the poloidal flux function (with $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}$ the poloidal flux), the toroidal angle variable, and $\unicode[STIX]{x1D717}$ the poloidal angle variable in straight field line coordinates ( $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}=I/qR^{2}$ ), and where $q=q(\unicode[STIX]{x1D713})$ is the safety factor and $I=I(\unicode[STIX]{x1D713})=RB_{t}$ . Here $R$ is the major radius and $B_{t}$ is the toroidal magnetic field. The Clebsch representation uses $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D701}-q\unicode[STIX]{x1D717}$ .

The Fokker–Planck equation written in the spatial and velocity space variables $\boldsymbol{r},\boldsymbol{v},t$ is

(3.2) $$\begin{eqnarray}\displaystyle {\dot{f}}=\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}t+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f+(Ze/M)(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}+c^{-1}\boldsymbol{v}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}f=C, & & \displaystyle\end{eqnarray}$$

where $f=f(\boldsymbol{r},\boldsymbol{v},t)$ is the distribution function, $\boldsymbol{E}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}$ is the electric field in the electrostatic limit, and $C$ denotes the Fokker–Planch collision operator. Also, $Z$ is the charge number and $M$ is the mass, with $e$ the charge on a proton and $c$ the speed of light. The mass ratio difference implies that it is usually convenient to think of the kinetic equation as being for the ions as electrons more closely follow field lines.

The axisymmetry of the magnetic field means that it will at times be convenient to employ the canonical angular momentum ( $-Ze\unicode[STIX]{x1D713}_{\ast }/c$ ) in the form

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{\ast }\equiv \unicode[STIX]{x1D713}-(Mc/Ze)R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}=\unicode[STIX]{x1D713}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-Iv_{\Vert }/\unicode[STIX]{x1D6FA}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}=ZeB/Mc$ and $v_{\Vert }=\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{v}$ . A departure of $\unicode[STIX]{x1D6F7}$ from axisymmetry means that $\unicode[STIX]{x1D713}_{\ast }$ is no longer a constant of the motion since

(3.4) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D713}}_{\ast }\equiv c\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701})\boldsymbol{\cdot }\boldsymbol{a}=0$ for an arbitrary vector $\boldsymbol{a}$ is used. The unperturbed vector potential $\text{}\underline{\boldsymbol{A}}$ is related to the flux function by $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}(\boldsymbol{r},t)=-R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\text{}\underline{\boldsymbol{A}}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D735}\times \text{}\underline{\boldsymbol{A}}$ .

Similarly, the total energy,

(3.5) $$\begin{eqnarray}\displaystyle E=Mv^{2}/2+Ze\unicode[STIX]{x1D6F7}, & & \displaystyle\end{eqnarray}$$

is not a constant of the motion because of the time variation of $\unicode[STIX]{x1D6F7}$ , which gives

(3.6) $$\begin{eqnarray}\displaystyle {\dot{E}}=Ze\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}/\unicode[STIX]{x2202}t. & & \displaystyle\end{eqnarray}$$

The Fokker–Planck kinetic equation is the starting point of all gyrokinetic (and drift kinetic) treatments.

4 Gyrokinetic orderings and assumptions

The gyrokinetic orderings assume a strongly magnetized plasma so the ion gyroradius $\unicode[STIX]{x1D70C}_{i}$ is small compared to the minor radius $a$ of the device, that is,

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{i}\ll a, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D70C}_{i}=v_{i}/\unicode[STIX]{x1D6FA}_{i}$ , $v_{i}=(2T_{i}/M_{i})^{1/2}$ the ion thermal speed and $\unicode[STIX]{x1D6FA}_{i}=ZeB/M_{i}c$ the ion cyclotron frequency. The minor radius $a$ is assumed to be a typical radial scale length, and is normally assumed much larger than the poloidal gyroradius $\unicode[STIX]{x1D70C}_{pi}=Bv_{i}/B_{p}\unicode[STIX]{x1D6FA}_{i}\simeq qR\unicode[STIX]{x1D70C}_{i}/a$ . The electrons are even more strongly magnetized.

In addition, for turbulent and collisional transport in tokamaks the frequencies, $\unicode[STIX]{x1D714}$ , of interest are normally assumed to be small compared to the ion cyclotron frequency. The typical frequencies of interest are near the diamagnetic drift frequency

(4.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D714}_{\ast }\sim k_{\bot }\unicode[STIX]{x1D70C}_{i}v_{i}/a & & \displaystyle\end{eqnarray}$$

associated with drift waves, where $k_{\bot }$ is a typical perpendicular wavenumber. The ion and electron drift frequencies are comparable as $\unicode[STIX]{x1D70C}_{i}v_{i}\sim \unicode[STIX]{x1D70C}_{e}v_{e}$ , with $\unicode[STIX]{x1D70C}_{e}$ and $v_{e}$ the electron gyroradius and thermal speed.

In the core of a tokamak collisions are normally weak so it is convenient to assume that the poloidal gyroradius is small compared to the ion mean free path $\unicode[STIX]{x1D706}=v_{i}/\unicode[STIX]{x1D708}_{ii}$ ,

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{pi}\ll \unicode[STIX]{x1D706}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{ii}$ is the ion–ion collision frequency, and ions and electrons have comparable mean free paths so no distinction is made.

Gyrokinetics allows

(4.4) $$\begin{eqnarray}\displaystyle k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1 & & \displaystyle\end{eqnarray}$$

by expanding in the small parameter

(4.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\sim \unicode[STIX]{x1D70C}_{i}/a\sim \unicode[STIX]{x1D714}_{\ast }/\unicode[STIX]{x1D6FA}_{i}\sim \unicode[STIX]{x1D70C}_{pi}/\unicode[STIX]{x1D706}. & & \displaystyle\end{eqnarray}$$

Moreover, $k_{\bot }\unicode[STIX]{x1D70C}_{e}\sim 1$ can be allowed so electron temperature gradient modes can be treated as well. In addition, the fluctuations are assumed to have parallel wavenumbers $k_{\Vert }$ such that

(4.6) $$\begin{eqnarray}\displaystyle k_{\Vert }qR\sim 1, & & \displaystyle\end{eqnarray}$$

where $qR$ is referred to as the connection length. Notice that $\unicode[STIX]{x1D70C}_{pi}/\unicode[STIX]{x1D706}=qR\unicode[STIX]{x1D70C}_{i}/\unicode[STIX]{x1D706}a$ implies the general collisional ordering $\unicode[STIX]{x1D706}\sim qR$ is allowed for both ions and electrons. These orderings are compatible with the banana regime ordering $\unicode[STIX]{x1D706}\gg qR$ in the core of a tokamak in which parallel streaming dominates over collisions. Moreover, they permit charges to stream faster along a field line than they magnetically drift poloidally in a flux surface, that is, $qR/v_{i}\ll a^{2}/\unicode[STIX]{x1D70C}_{i}v_{i}$ , since $\unicode[STIX]{x1D70C}_{pi}/a\ll 1$ . Only on the longer collisional or turbulent transport scale do charges transport across flux surfaces.

It is normally useful to write the electrostatic potential as

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}=\bar{\unicode[STIX]{x1D6F7}}(\boldsymbol{r},t)+\tilde{\unicode[STIX]{x1D6F7}}(\boldsymbol{r},t), & & \displaystyle\end{eqnarray}$$

where $\bar{\unicode[STIX]{x1D6F7}}$ and $\tilde{\unicode[STIX]{x1D6F7}}$ , respectively, are slowly and rapidly varying functions of space and time. For a sonic ( ${\sim}v_{i}$ ) electric field drift $Ze\bar{\unicode[STIX]{x1D6F7}}/T\sim 1/\unicode[STIX]{x1D6FF}\gg 1$ . Normally, however, the drift due to the electric field is at the level of the diamagnetic drift speed ( ${\sim}\unicode[STIX]{x1D6FF}v_{i}$ ) giving $Ze\bar{\unicode[STIX]{x1D6F7}}/T\sim 1$ .

It is also often convenient to employ scale separation to split the distribution function up into slow and fast spatio-temporal portions by writing

(4.8) $$\begin{eqnarray}\displaystyle f=\bar{f}(\boldsymbol{r},\boldsymbol{v},t)+\tilde{f}(\boldsymbol{r},\boldsymbol{v},t). & & \displaystyle\end{eqnarray}$$

Typically $\unicode[STIX]{x2202}\tilde{f}/\unicode[STIX]{x2202}t\sim \unicode[STIX]{x1D714}_{\ast }\tilde{f}$ and $\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t\sim \unicode[STIX]{x1D714}_{\ast }\tilde{\unicode[STIX]{x1D6F7}}$ , with $\bar{f}$ and $\bar{\unicode[STIX]{x1D6F7}}$ varying approximately on the transport scale of gyroDohm diffusion $D_{gB}\sim \unicode[STIX]{x1D70C}_{i}^{2}v_{i}/a$ , that is, $\bar{f}^{-1}\unicode[STIX]{x2202}\bar{f}/\unicode[STIX]{x2202}t\sim D_{gB}/a^{2}$ .

Moreover, it is useful to order

(4.9) $$\begin{eqnarray}\displaystyle \tilde{f}/\bar{f}\sim Ze\tilde{\unicode[STIX]{x1D6F7}}/T\sim 1/k_{\bot }a\ll 1, & & \displaystyle\end{eqnarray}$$

which is consistent with the so-called mixing length estimate $\unicode[STIX]{x1D735}_{\bot }\,\bar{f}\sim \unicode[STIX]{x1D735}_{\bot }\,\tilde{f}$ . The preceding ordering for $\tilde{f}/\bar{f}\sim Ze\tilde{\unicode[STIX]{x1D6F7}}/T$ was introduced by Dimits, Lodestro & Dubin (Reference Dimits, Lodestro and Dubin1992). Using it gives the fluctuating electric drift to be of order

(4.10) $$\begin{eqnarray}\displaystyle cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}\sim \unicode[STIX]{x1D6FF}v_{i}. & & \displaystyle\end{eqnarray}$$

Also, rapid temporal variation on the diamagnetic drift frequency scale, and slow temporal variation on the gyroBohm transport scale gives $(\unicode[STIX]{x2202}\bar{f}/\unicode[STIX]{x2202}t)/(\unicode[STIX]{x2202}\tilde{f}/\unicode[STIX]{x2202}t)\sim (\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t)/(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t)\sim \unicode[STIX]{x1D6FF}$ for

(4.11) $$\begin{eqnarray}\displaystyle Ze\bar{\unicode[STIX]{x1D6F7}}/T\sim 1. & & \displaystyle\end{eqnarray}$$

5 Gyrokinetic variables

When transforming to gyrokinetic variables there are a few choices for the velocity space variables, but the usual choice for the spatial variation is the mixed variable that relates the guiding centre or gyrocentre location $\boldsymbol{R}$ to the particle position $\boldsymbol{r}$ , namely,

(5.1) $$\begin{eqnarray}\displaystyle \boldsymbol{R}=\boldsymbol{r}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}=\boldsymbol{r}-\unicode[STIX]{x1D746}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D746}$ the gyroradius vector. The velocity variable is defined in cylindrical velocity space variables $v_{\bot },\unicode[STIX]{x1D711},v_{\Vert }$ such that the charge spirals about the magnetic field line as it streams along it to lowest order. Then

(5.2) $$\begin{eqnarray}\displaystyle \boldsymbol{v}=v_{\bot }(\boldsymbol{e}_{\unicode[STIX]{x1D713}}\cos \unicode[STIX]{x1D711}+\boldsymbol{e}_{\times }\sin \unicode[STIX]{x1D711})+v_{\Vert }\boldsymbol{b}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D711}$ the gyrophase, and $\boldsymbol{e}_{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}/|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|$ and $\boldsymbol{e}_{\times }=\boldsymbol{b}\times \boldsymbol{e}_{\unicode[STIX]{x1D713}}$ orthonormal unit vectors.

There are three main reasons why this gyrokinetic change of variables is useful. The first is that to lowest order the guiding centre tries to stream along the magnetic field, that is,

(5.3) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{R}}\simeq \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{r}+(Ze/Mc)(\boldsymbol{v}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}(\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b})=\boldsymbol{v}-\boldsymbol{v}_{\bot }=v_{\Vert }\boldsymbol{b}. & & \displaystyle\end{eqnarray}$$

To obtain the preceding, the velocity space variable result $\unicode[STIX]{x1D735}_{v}\boldsymbol{v}=\overset{\leftrightarrow }{I}$ is employed, just as in spatial variables $\unicode[STIX]{x1D735}\boldsymbol{r}=\overset{\leftrightarrow }{I}$ . A convenient form of the unit dyad is $\overset{\leftrightarrow }{I}=\boldsymbol{e}_{\unicode[STIX]{x1D713}}\boldsymbol{e}_{\unicode[STIX]{x1D713}}+\boldsymbol{e}_{\times }\boldsymbol{e}_{\times }+\boldsymbol{b}\boldsymbol{b}$ .

The preceding means that the gyrokinetic variable $\boldsymbol{R}$ allows the ordering

(5.4) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}f\sim \unicode[STIX]{x1D6FA}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}f. & & \displaystyle\end{eqnarray}$$

Consequently, unlike drift kinetics (Hazeltine Reference Hazeltine1973), gyrokinetics is able to retain short perpendicular wavelengths by allowing $\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{v}_{\bot }\sim \unicode[STIX]{x1D6FA}$ .

It is convenient to employ scale separation to assume $\boldsymbol{B}$ , $\bar{\unicode[STIX]{x1D6F7}}$ and $\bar{f}$ are slow functions of space such that Taylor expansions may be employed to obtain

(5.5) $$\begin{eqnarray}\displaystyle \boldsymbol{B}(\boldsymbol{r})=\boldsymbol{B}(\boldsymbol{R})+\unicode[STIX]{x1D746}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{B}+\cdots \simeq \boldsymbol{B}(\boldsymbol{R}), & & \displaystyle\end{eqnarray}$$

(5.6) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D6F7}}(\boldsymbol{r})\simeq \bar{\unicode[STIX]{x1D6F7}}(\boldsymbol{R})+\unicode[STIX]{x1D746}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\unicode[STIX]{x1D6F7}}+\cdots \simeq \bar{\unicode[STIX]{x1D6F7}}(\boldsymbol{R}), & & \displaystyle\end{eqnarray}$$

and

(5.7) $$\begin{eqnarray}\displaystyle \bar{f}(\boldsymbol{R})\simeq \bar{f}(\boldsymbol{r})-\unicode[STIX]{x1D746}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{f}+\cdots \,, & & \displaystyle\end{eqnarray}$$

where the leading correction to the expansion of $\bar{f}$ is needed to treat diamagnetic effects.

The fluctuating quantities $\widetilde{\unicode[STIX]{x1D6F7}}$ and $\widetilde{f}$ contain rapid spatial variation and need not be Taylor expanded. This feature is another important difference between gyrokinetics and drift kinetics. For example, if

(5.8) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D6F7}}\propto \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}=\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{R}+\unicode[STIX]{x1D746})}, & & \displaystyle\end{eqnarray}$$

then a gyroaverage at fixed gyrocenter location $\boldsymbol{R}$ ,

(5.9) $$\begin{eqnarray}\displaystyle \langle \cdots \rangle _{R}=(2\unicode[STIX]{x03C0})^{-1}\oint \text{d}\unicode[STIX]{x1D711}(...), & & \displaystyle\end{eqnarray}$$

may be employed to obtain

(5.10) $$\begin{eqnarray}\displaystyle \langle \exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r})\rangle _{R}=\langle \exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746})\rangle _{R}\exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R})=J_{0}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA})\exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}). & & \displaystyle\end{eqnarray}$$

As a result, the gyroaverage of the electrostatic potential holding the guiding centre fixed differs from the electrostatic potential at fixed particle location, that is, $\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}\neq \tilde{\unicode[STIX]{x1D6F7}}(\boldsymbol{r})$ .

The third key reason that the gyrokinetic change of variables is so useful is found by working a bit harder. Assuming a steady state or weakly time varying magnetic field,

(5.11) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{R}}=v_{\Vert }\boldsymbol{b}-\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b})\times \boldsymbol{v}-cB^{-1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}\times \boldsymbol{b}. & & \displaystyle\end{eqnarray}$$

Using $\langle \boldsymbol{v}\rangle _{R}=v_{\Vert }\boldsymbol{b}$ and $\langle \boldsymbol{v}\boldsymbol{v}\rangle _{R}=(v_{\bot }^{2}/2)(\overset{\leftrightarrow }{I}-\boldsymbol{b}\boldsymbol{b})+v_{\Vert }^{2}\boldsymbol{b}\boldsymbol{b}$ gives the important result

(5.12) $$\begin{eqnarray}\displaystyle \langle \dot{\boldsymbol{R}}\rangle _{R}=v_{\Vert }\boldsymbol{b}-\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b})\times \boldsymbol{v}\rangle _{R}+cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}=(v_{\Vert }+u_{\Vert })\boldsymbol{b}+\boldsymbol{v}_{M}+\langle \boldsymbol{v}_{E}\rangle _{R}, & & \displaystyle\end{eqnarray}$$

where the gyroaveraged electric drift is

(5.13) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{E}\rangle _{R}=cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle \unicode[STIX]{x1D6F7}\rangle _{R} & & \displaystyle\end{eqnarray}$$

and the magnetic drift is as in drift kinetics

(5.14) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{M}=\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b}\times [v_{\Vert }^{2}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}+(v_{\bot }^{2}/2)\unicode[STIX]{x1D735}\ell nB], & & \displaystyle\end{eqnarray}$$

while the parallel velocity correction is

(5.15) $$\begin{eqnarray}\displaystyle u_{\Vert }=-(v_{\bot }^{2}/2\unicode[STIX]{x1D6FA})\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}. & & \displaystyle\end{eqnarray}$$

To more easily obtain $u_{\Vert }$ and the $\unicode[STIX]{x1D735}B$ and curvature drifts, it is useful to notice that

(5.16) $$\begin{eqnarray}\displaystyle -\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b})\times \boldsymbol{v}\rangle _{R}=(v_{\bot }^{2}/2)\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b})+\unicode[STIX]{x1D6FA}^{-1}(v_{\Vert }^{2}-v_{\bot }^{2}/2)\boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}), & & \displaystyle\end{eqnarray}$$

and then use $\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b})=\unicode[STIX]{x1D6FA}^{-1}(\unicode[STIX]{x1D735}\times \boldsymbol{b}+\boldsymbol{b}\times \unicode[STIX]{x1D735}\ell nB)$ and $(\overset{\leftrightarrow }{I}-\boldsymbol{b}\boldsymbol{b})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}=\boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b})$ to find

(5.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b})=\unicode[STIX]{x1D6FA}^{-1}[\boldsymbol{b}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}+\boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}+\unicode[STIX]{x1D735}\ell nB)]. & & \displaystyle\end{eqnarray}$$

The combination

(5.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D713}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713} & & \displaystyle\end{eqnarray}$$

appearing in the canonical angular momentum is the radial gyrokinetic variable, with the poloidal and toroidal gyrokinetic variables similarly defined as

(5.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D717}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}, & & \displaystyle\end{eqnarray}$$

and

(5.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D701}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}. & & \displaystyle\end{eqnarray}$$

The drift kinetic canonical angular momentum is

(5.21) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D713}}_{\ast }=\unicode[STIX]{x1D713}-Iv_{\Vert }/\unicode[STIX]{x1D6FA}, & & \displaystyle\end{eqnarray}$$

a constant of the motion of the drift kinetic equation. Occasionally, it is convenient to use the canonical angular momentum $\unicode[STIX]{x1D713}_{\ast }$ as the radial variable rather than the gyrokinetic variable $\unicode[STIX]{x1D6F9}$ (Kagan & Catto Reference Kagan and Catto2008, Reference Kagan and Catto2010).

Analytically, it is usually most convenient to use the velocity space variables of total energy,

(5.22) $$\begin{eqnarray}\displaystyle E=Mv^{2}/2+Ze\unicode[STIX]{x1D6F7}, & & \displaystyle\end{eqnarray}$$

magnetic moment,

(5.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}=Mv_{\bot }^{2}/2B, & & \displaystyle\end{eqnarray}$$

and gyrophase $\unicode[STIX]{x1D711}$ .

In these variables a so-called full f form of the gyrokinetic equation can be obtained by gyroaveraging

(5.24) $$\begin{eqnarray}\displaystyle {\dot{f}}=\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}t+\dot{\boldsymbol{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}f+{\dot{E}}\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}E+\dot{\unicode[STIX]{x1D707}}\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}\unicode[STIX]{x1D707}+\dot{\unicode[STIX]{x1D711}}\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}. & & \displaystyle\end{eqnarray}$$

For slow temporal evolution and global variation of the magnetic field

(5.25) $$\begin{eqnarray}\displaystyle & & \displaystyle {\dot{E}}=Ze\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}/\unicode[STIX]{x2202}t,\end{eqnarray}$$

(5.26) $$\begin{eqnarray}\displaystyle & & \displaystyle \dot{\unicode[STIX]{x1D707}}=-\unicode[STIX]{x1D707}B^{-1}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B-Mv_{\Vert }B^{-1}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{v}-ZeB^{-1}\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7},\end{eqnarray}$$

and, differentiating $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{e}_{\unicode[STIX]{x1D713}}=v_{\bot }\cos \unicode[STIX]{x1D711}$ and $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{e}_{\times }=v_{\bot }\sin \unicode[STIX]{x1D711}$ ,

(5.27) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D711}}=-\unicode[STIX]{x1D6FA}-\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{e}_{\unicode[STIX]{x1D713}}\boldsymbol{\cdot }\boldsymbol{e}_{\times }-v_{\Vert }v_{\bot }^{-2}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}\times \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{v}+(Ze/Mv_{\bot }^{2})\boldsymbol{v}\;\boldsymbol{\cdot }\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}. & & \displaystyle\end{eqnarray}$$

These last three expressions are the same as for drift kinetics. However, gyroaveraging holding $\boldsymbol{R}$ fixed instead of $\boldsymbol{r}$ , gives the lowest-order expressions

(5.28) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle {\dot{E}}\rangle _{R}=Ze\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D6F7}\rangle _{R}/\unicode[STIX]{x2202}t,\end{eqnarray}$$

(5.29) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle \dot{\unicode[STIX]{x1D707}}\rangle _{R}=0,\end{eqnarray}$$

and

(5.30) $$\begin{eqnarray}\displaystyle \langle \dot{\unicode[STIX]{x1D711}}\rangle _{R}=-\unicode[STIX]{x1D6FA}, & & \displaystyle\end{eqnarray}$$

where the lowest-order result $\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}=(\unicode[STIX]{x2202}\boldsymbol{v}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711})\boldsymbol{\cdot }\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}|_{R}$ and $Ze\bar{\unicode[STIX]{x1D6F7}}/T\sim 1$ are used, and $\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{e}_{\unicode[STIX]{x1D713}}\boldsymbol{\cdot }\boldsymbol{e}_{\times }+v_{\Vert }v_{\bot }^{-2}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}\times \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{v}\rangle _{R}\simeq v_{\Vert }(2\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{e}_{\unicode[STIX]{x1D713}}\boldsymbol{\cdot }\boldsymbol{e}_{\times }+\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b})/2\sim v_{\Vert }/R\ll \unicode[STIX]{x1D6FA}$ varies on the scale of the magnetic field so is small. Expanding f in $\unicode[STIX]{x1D6FF}$ by writing $f=f^{(0)}+f^{(1)}+\cdots \,,$ taking $-\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}f^{(0)}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$ to lowest order, and then gyroaveraging the next-order equation to annihilate the $f^{(1)}$ term, gives the simplest full f gyrokinetic equation to be

(5.31) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\unicode[STIX]{x2202}f^{(0)}/\unicode[STIX]{x2202}t+(v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{M}+\langle \boldsymbol{v}_{E}\rangle _{R})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}f^{(0)}+Ze(\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D6F7}\rangle _{R}/\unicode[STIX]{x2202}t)\unicode[STIX]{x2202}f^{(0)}/\unicode[STIX]{x2202}E=\langle C\{\,f^{(0)}\}\rangle _{R}\\ \quad \unicode[STIX]{x1D714}_{\ast }\quad ~v_{i}/qR\quad k_{\bot }\unicode[STIX]{x1D70C}_{i}v_{i}/R\quad ~\unicode[STIX]{x1D6FF}v_{i}/a\qquad \qquad \quad \quad \unicode[STIX]{x1D714}_{\ast }/k_{\bot }a\quad \qquad \qquad \quad \unicode[STIX]{x1D708}.\end{array}\right\}\quad & & \displaystyle\end{eqnarray}$$

Unfortunately, this direct approach is risky since the result contains terms of a few different orders as noted below the equation (recall $k_{\bot }a\gg 1$ and $\unicode[STIX]{x1D708}\ll \unicode[STIX]{x1D714}_{\ast }$ ), and it ignores higher-order corrections from $f^{(1)}$ and the gyrokinetic change of variables. In addition, most of $f^{(0)}$ is normally a Maxwellian. Consequently, the derivation of the gyrokinetic equation will be refined in the next sections by being more clever and deriving what is referred to as a local or $\unicode[STIX]{x1D6FF}f=\tilde{f}$ gyrokinetic equation. The derivation takes advantage of the scale separation between the fluctuations and the background.

6 Gyrokinetics for the slowly varying portion of the distribution function

When the unperturbed distribution function is slowly varying in time and space with $\unicode[STIX]{x2202}\bar{f}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$ , only the lowest-order equation

(6.1) $$\begin{eqnarray}\displaystyle \langle \dot{\bar{f}}\rangle _{R}\simeq v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}=C\{\langle \bar{f}\rangle _{r}\}\simeq C\{\,\bar{f}\}, & & \displaystyle\end{eqnarray}$$

need be solved, as long as $\unicode[STIX]{x1D708}\bar{f}\gg \unicode[STIX]{x2202}\bar{f}/\unicode[STIX]{x2202}t$ , $\langle \boldsymbol{v}_{E}\rangle _{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}/v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}\sim \boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}/v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}\sim \unicode[STIX]{x1D70C}_{pi}/a\ll 1$ and

(6.2) $$\begin{eqnarray}\displaystyle \bar{f}(\boldsymbol{R},\bar{E},\unicode[STIX]{x1D707},t)\simeq \bar{f}(\boldsymbol{r},\bar{E},\unicode[STIX]{x1D707},t)+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f+\cdots \,. & & \displaystyle\end{eqnarray}$$

The preceding implies that $v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\,\bar{f}\simeq v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}$ and $\langle C\{\,\bar{f}\}\rangle _{R}\simeq \langle C\{\,\bar{f}\}\rangle _{r}=C\{\langle \,\bar{f}\rangle _{r}\}\simeq C\{\,\bar{f}\}$ . A subscript r on a gyroaverage indicates that it is to be performed at fixed $\boldsymbol{r}$ rather than fixed $\boldsymbol{R}$ . As $\tilde{\unicode[STIX]{x1D6F7}}$ is small compared to $\bar{\unicode[STIX]{x1D6F7}}$ it will be convenient at times to use the alternate energy variable

(6.3) $$\begin{eqnarray}\displaystyle \bar{E}=Mv^{2}/2+Ze\bar{\unicode[STIX]{x1D6F7}}. & & \displaystyle\end{eqnarray}$$

Using the lowest-order form of the kinetic equation $v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}=C\{\,\bar{f}\}$ , either a flux surface average,

(6.4) $$\begin{eqnarray}\displaystyle \langle \cdots \rangle =\left[\left.\oint \text{d}\unicode[STIX]{x1D717}\text{d}\unicode[STIX]{x1D701}(...)/\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\right]\right/\left[\oint \text{d}\unicode[STIX]{x1D717}\text{d}\unicode[STIX]{x1D701}/\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\right], & & \displaystyle\end{eqnarray}$$

for the passing, or a transit average,

(6.5) $$\begin{eqnarray}\displaystyle \overline{(...)}=\left[\left.\oint \text{d}\unicode[STIX]{x1D717}\text{d}\unicode[STIX]{x1D701}(...)/v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\right]\right/\left[\oint \text{d}\unicode[STIX]{x1D717}\text{d}\unicode[STIX]{x1D701}/v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\right], & & \displaystyle\end{eqnarray}$$

for the trapped, can be employed to annihilate the parallel streaming term. It is convenient to retain the $\unicode[STIX]{x1D701}$ integrals in the transit average, and also notice $\overline{v_{\Vert }/B}=1/\langle B/v_{\Vert }\rangle$ for the passing, with $\overline{v_{\Vert }/B}=0$ for the trapped. The constraints obtained for the passing by flux surface averaging,

(6.6) $$\begin{eqnarray}\displaystyle \langle Bv_{\Vert }^{-1}C\{\,\bar{f}\}\rangle =0, & & \displaystyle\end{eqnarray}$$

and trapped by transit averaging,

(6.7) $$\begin{eqnarray}\displaystyle \overline{C\{\,\bar{f}\}}=0, & & \displaystyle\end{eqnarray}$$

must be satisfied by the lowest-order distribution function as is shown next.

To verify that these can only be satisfied by a Maxwellian, multiply the lowest-order kinetic equation by $\ell n\bar{f}$ and integrate over all velocity space and flux surface average to obtain the Boltzmann H theorem on a flux surface,

(6.8) $$\begin{eqnarray}\displaystyle \left\langle \int \text{d}^{3}v\ell n\bar{f}C\{\,\bar{f}\}\right\rangle =0. & & \displaystyle\end{eqnarray}$$

The preceding form for the H theorem follows from $\text{d}^{3}v\rightarrow 2\unicode[STIX]{x03C0}\text{d}\bar{E}\text{d}\unicode[STIX]{x1D707}B/v_{\Vert }$ and

(6.9) $$\begin{eqnarray}\displaystyle \left\langle \int \text{d}^{3}vv_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\bar{f}\ell n\bar{f}-\bar{f})\right\rangle =\left\langle \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left[B^{-1}\int \text{d}^{3}vv_{\Vert }(\bar{f}\ell n\bar{f}-\bar{f})\right]\right\rangle =0 & & \displaystyle\end{eqnarray}$$

because the flux surface average of the divergence of any vector $\boldsymbol{Q}$ ,

(6.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Q}=\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\left[\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left(\displaystyle \frac{\boldsymbol{Q}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}}\right)+\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}\left(\frac{\boldsymbol{Q}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}}{\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\left(\displaystyle \frac{\boldsymbol{Q}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}}{\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}}\right)\right], & & \displaystyle\end{eqnarray}$$

gives

(6.11) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Q}\rangle =\displaystyle \frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}(V^{\prime }\langle \boldsymbol{Q}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle ), & & \displaystyle\end{eqnarray}$$

with $V^{\prime }=\oint \text{d}\unicode[STIX]{x1D717}\text{d}\unicode[STIX]{x1D701}/\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ and where $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ is independent of $\unicode[STIX]{x1D701}$ . The entropy production only vanishes for a local Maxwellian $f_{M}$ , for which $C\{\,\bar{f}\}=0$ and, therefore, $v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{f}=0$ . Consequently, the steady state solution for a stationary confined plasma with well-defined flux or total pressure surfaces is

(6.12) $$\begin{eqnarray}\displaystyle \bar{f}=f_{M}(\unicode[STIX]{x1D713},\bar{E})=f_{M}=\unicode[STIX]{x1D702}(\unicode[STIX]{x1D713})[M/2\unicode[STIX]{x03C0}T(\unicode[STIX]{x1D713})]^{3/2}\text{e}^{-\bar{E}/T(\unicode[STIX]{x1D713})}=n(M/2\unicode[STIX]{x03C0}T)^{3/2}\text{e}^{-Mv^{2}/2T}, & & \displaystyle\end{eqnarray}$$

with the density $n$ and pseudo-density $\unicode[STIX]{x1D702}$ related by

(6.13) $$\begin{eqnarray}\displaystyle n=\unicode[STIX]{x1D702}(\unicode[STIX]{x1D713})\text{e}^{-Ze\bar{\unicode[STIX]{x1D6F7}}(\unicode[STIX]{x1D713})/T(\unicode[STIX]{x1D713})}. & & \displaystyle\end{eqnarray}$$

In the core of a tokamak $n$ , $T$ and $\bar{\unicode[STIX]{x1D6F7}}$ are lowest-order flux functions. The collision operator must be Galilean invariant so a drifting Maxwellian is also allowed. However, without large momentum input, the ion drift velocity is small and at the diamagnetic level ${\sim}\unicode[STIX]{x1D6FF}v_{i}$ so a stationary Maxwellian can be employed.

The preceding demonstration is valid for ions without any modification since $C$ is then the ion–ion collision operator. For electrons both like and unlike electron–ion collisions must be retained. However, momentum relaxation occurs on the electron–ion time scale so friction rapidly relaxes the electrons to the same drift velocity as the ions. A speed expansion of the electron–ion collision operator can then be performed and the Galilean invariance of the like collision operator used to determine that the electrons are also Maxwellian to lowest order.

The next-order correction to the Maxwellian in a tokamak is found by using axisymmetry. To find this correction it is convenient to define

(6.14) $$\begin{eqnarray}\displaystyle \bar{f}_{\ast }=\unicode[STIX]{x1D702}(\unicode[STIX]{x1D713}_{\ast })[M/2\unicode[STIX]{x03C0}T(\unicode[STIX]{x1D713}_{\ast })]^{3/2}\text{e}^{-\bar{E}/T(\unicode[STIX]{x1D713}_{\ast })}, & & \displaystyle\end{eqnarray}$$

so that Taylor expansion gives

(6.15) $$\begin{eqnarray}\displaystyle \overline{f_{\ast }}(\unicode[STIX]{x1D713}_{\ast },\bar{E})\simeq f_{M}(\unicode[STIX]{x1D713}\bar{E})-(Mc/Ze)R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}+\cdots \,. & & \displaystyle\end{eqnarray}$$

For an axisymmetric steady state

(6.16) $$\begin{eqnarray}\displaystyle \dot{\overline{f}}_{\ast }(\unicode[STIX]{x1D713}_{\ast },\bar{E})=0, & & \displaystyle\end{eqnarray}$$

since $\dot{\unicode[STIX]{x1D713}}_{\ast }=0$ and $\dot{\bar{E}}=0$ . Consequently, letting

(6.17) $$\begin{eqnarray}\displaystyle \overline{f}=\overline{f}_{\ast }(\unicode[STIX]{x1D713}_{\ast },\bar{E})+\overline{h}, & & \displaystyle\end{eqnarray}$$

and using $\langle C_{1}\{(\unicode[STIX]{x1D713}_{\ast }-\unicode[STIX]{x1D713})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}\rangle _{R}\simeq -C_{1}\{(Iv_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}$ , the gyrokinetic equation for $\overline{h}$ is

(6.18) $$\begin{eqnarray}\displaystyle \langle \dot{\bar{f}}\rangle _{R}\simeq v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{h}=C_{1}\{\overline{h}-(Iv_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}, & & \displaystyle\end{eqnarray}$$

with $C_{1}$ denoting the linearized collision operator.

To lowest-order collisions are normally weak in the core, giving $v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{h}=0$ . Transit averaging the next-order equation for the trapped (subscript $t$ ) over a full bounce annihilates the streaming term leaving

(6.19) $$\begin{eqnarray}\displaystyle 0=\overline{C_{1}\{\overline{h}_{t}-(Iv_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}}=C_{1}\{\overline{h}_{t}-(I\overline{v_{\Vert }}/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}=C_{1}\{\overline{h}_{t}\}. & & \displaystyle\end{eqnarray}$$

As there is no drive,

(6.20) $$\begin{eqnarray}\displaystyle \overline{h}_{t}=0. & & \displaystyle\end{eqnarray}$$

A flux surface average is used to annihilate the streaming term in the next-order equation for the passing (subscript $p$ ) leaving

(6.21) $$\begin{eqnarray}\displaystyle \langle Bv_{\Vert }^{-1}C_{1}\{\bar{h}_{p}-(Iv_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}\rangle =0. & & \displaystyle\end{eqnarray}$$

The solution of this constraint is the neoclassical response $\bar{h}_{p}=\bar{h}_{p}(\unicode[STIX]{x1D713},\bar{E},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})$ responsible for banana regime transport. It must be odd in $v_{\Vert }=\unicode[STIX]{x1D70E}|v_{\Vert }|$ like the drive. As it depends only on $\unicode[STIX]{x1D713}$ , $\bar{E}$ , and $\unicode[STIX]{x1D707}$ it will be shown later that it does not matter in the gyrokinetic equation for the fluctuations even though $\bar{h}_{p}\sim (Iv_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\sim (Mc/Ze)R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$ .

Before going on to derive the gyrokinetic equation for the fluctuations some particle transport and ambipolarity considerations are worth stressing.

7 Ambipolarity and particle transport considerations

The convenient way to evaluate the collisional and turbulent particle fluxes is to use the $R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ component of the momentum conservation equation

(7.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}(Mn\boldsymbol{V})/\unicode[STIX]{x2202}t+Zen(\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}-c^{-1}\boldsymbol{V}\times \boldsymbol{B})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(M\int \text{d}^{3}v\boldsymbol{v}\boldsymbol{v}f)=M\int \text{d}^{3}v\boldsymbol{v}C. & & \displaystyle\end{eqnarray}$$

Forming the $R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ component and including coarse grain averages in time and radius with the flux surface average gives

(7.2) $$\begin{eqnarray}\displaystyle \displaystyle \frac{Ze}{c}\langle n\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg} & = & \displaystyle Ze\left\langle n\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\right\rangle _{cg}+\displaystyle \frac{1}{V^{\prime }}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }\langle MR^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}vf\boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg}\nonumber\\ \displaystyle & & \displaystyle \quad -\,M\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}v\boldsymbol{v}C\rangle _{cg},\end{eqnarray}$$

where the density and flow velocity are defined as $n=\int \text{d}^{3}vf$ and $n\boldsymbol{V}=\int \text{d}^{3}v\boldsymbol{v}f$ , and

(7.3) $$\begin{eqnarray}\displaystyle \langle \cdots \rangle _{cg}=\displaystyle \frac{1}{2\unicode[STIX]{x1D70F}}\int _{t-\unicode[STIX]{x1D70F}}^{t+\unicode[STIX]{x1D70F}}\text{d}t\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\int _{\unicode[STIX]{x1D713}-\unicode[STIX]{x1D6FE}}^{\unicode[STIX]{x1D713}+\unicode[STIX]{x1D6FE}}\text{d}\unicode[STIX]{x1D713}\langle \cdots \rangle , & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D714}_{\ast }\unicode[STIX]{x1D70F}\gg 1$ and $\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{e}_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D6FE}\gg |\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|$ . The coarse grain average removes very fast time and very short radial scale variation. The species pressure $p=(M/3)\int \text{d}^{3}vv^{2}f=nT$ does not contribute to $\langle MR^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}vf\boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle$ , making the term small. Indeed, the off-diagonal stress tensor and time derivative only enter at higher order when the radial flux of toroidal angular momentum is evaluated. Therefore, only the turbulent and collisional fluxes survive in

(7.4) $$\begin{eqnarray}\displaystyle Ze\langle n\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg}=Zec\left\langle n\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\right\rangle _{cg}-Mc\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}v\boldsymbol{v}C\rangle _{cg}. & & \displaystyle\end{eqnarray}$$

In a quasineutral plasma ( $\unicode[STIX]{x1D6F4}Zen=0$ ) this form makes it clear that turbulent and collisional particle transport are separately intrinsically ambipolar, and when added yield

(7.5) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{J}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle =0, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{J}$ is the current density and $\unicode[STIX]{x1D6F4}$ is the species sum. This ambipolarity condition is required in a quasineutral plasma for which $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}=0$ must satisfy the constraint $\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}\rangle =0$ with the radial current vanishing at the magnetic axis. Overall momentum conservation for Coulomb collisions requires

(7.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F4}\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }M\int \text{d}^{3}v\boldsymbol{v}C\rangle _{cg}=0, & & \displaystyle\end{eqnarray}$$

as first noted by Kovrizhnykh (Reference Kovrizhnykh1969) and Rutherford (Reference Rutherford1970). The form of $\langle n\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle$ implies that the turbulent particle flux is

(7.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{\text{turb}}=Zec\langle n\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}\rangle _{cg} & & \displaystyle\end{eqnarray}$$

and the collisional particle flux is

(7.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{\text{coll}}=-cM\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}v(\boldsymbol{v}_{\bot }+v_{\Vert }\boldsymbol{b})C\rangle _{cg}, & & \displaystyle\end{eqnarray}$$

with the $\boldsymbol{v}_{\bot }$ the contribution from classical collisions due to gyromotion and the $v_{\Vert }\boldsymbol{b}$ term the neoclassical portion caused by finite orbit departures from flux surfaces. The separate intrinsic ambipolarity of turbulent ( $\unicode[STIX]{x1D6F4}\unicode[STIX]{x1D6E4}_{\text{turb}}=0$ , using quasineutrality) and collisional ( $\unicode[STIX]{x1D6F4}\unicode[STIX]{x1D6E4}_{\text{coll}}=0$ , using collisional momentum conservation) particle transport in a tokamak was first demonstrated by Sugama et al. (Reference Sugama, Okamoto, Horton and Wakatani1996). Their treatment is fully electromagnetic, but some details relied on using a quasilinear treatment and Krook collision operator. These limitations where removed electrostatically by Parra & Catto (Reference Parra and Catto2009). Appendix A of Sugama & Horton (Reference Sugama and Horton1998) also showed intrinsic ambipolarity electromagnetically in a tokamak with sonic toroidal flow within a quasilinear treatment of the Onsager symmetry.

In a completely axisymmetric tokamak it is necessary to satisfy conservation of total toroidal angular momentum to higher order. Quasineutrality, ambipolarity and total collisional momentum conservation require that the lowest-order flux function behaviour of $\bar{\unicode[STIX]{x1D6F7}}$ adjusts to satisfy conservation of toroidal angular momentum

(7.9) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D6F4}M\langle R^{2}n\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\rangle +\displaystyle \frac{c}{V^{\prime }}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }\unicode[STIX]{x1D6F4}M\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}vf\boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg}=0, & & \displaystyle\end{eqnarray}$$

where the time derivative term is inserted for completeness. In the absence of a momentum source (as assumed here) the steady state equation reduces to the condition that the radial flux of toroidal angular momentum vanish on each flux surface,

(7.10) $$\begin{eqnarray}\displaystyle M\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}vf\boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg}=0, & & \displaystyle\end{eqnarray}$$

as there is no flux of radial momentum at the magnetic axis. The species sum is dropped as the electron contribution is usually negligible.

8 Gyrokinetic equation for the fluctuations

The conventional $\unicode[STIX]{x1D6FF}f=\tilde{f}$ gyrokinetic equation is really not an equation for $\unicode[STIX]{x1D6FF}f=\tilde{f}$ at all. Instead it is the gyrokinetic equation for the non-adiabatic portion of $\unicode[STIX]{x1D6FF}f=\tilde{f}$ obtained by letting

(8.1) $$\begin{eqnarray}\displaystyle f=f_{\ast }(\unicode[STIX]{x1D713}_{\ast },E)+h=\bar{f}+\tilde{f}, & & \displaystyle\end{eqnarray}$$

with

(8.2) $$\begin{eqnarray}\displaystyle h=\bar{h}+\tilde{h}, & & \displaystyle\end{eqnarray}$$

and $\tilde{h}$ the non-adiabatic response. Expanding the $E$ dependence about

(8.3) $$\begin{eqnarray}\displaystyle \bar{E}=Mv^{2}/2+Ze\bar{\unicode[STIX]{x1D6F7}}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717},t) & & \displaystyle\end{eqnarray}$$

for small $\tilde{\unicode[STIX]{x1D6F7}}$ , and expanding $\unicode[STIX]{x1D713}_{\ast }$ about $\unicode[STIX]{x1D713}$ yields

(8.4) $$\begin{eqnarray}\displaystyle f_{\ast }(\unicode[STIX]{x1D713}_{\ast },E)\simeq f_{M}(\unicode[STIX]{x1D713},\bar{E})-(Ze\tilde{\unicode[STIX]{x1D6F7}}/T)f_{M}-(Mc/Ze)R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}+\cdots \,. & & \displaystyle\end{eqnarray}$$

Then

(8.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f=\tilde{f}=\tilde{h}-(Ze\tilde{\unicode[STIX]{x1D6F7}}/T)f_{M}, & & \displaystyle\end{eqnarray}$$

with $-(Ze\tilde{\unicode[STIX]{x1D6F7}}/T)f_{M}$ the adiabatic part or Maxwell–Boltzmann response of $\tilde{f}$ , and

(8.6) $$\begin{eqnarray}\displaystyle \overline{f}=\overline{f_{\ast }}(\unicode[STIX]{x1D713}_{\ast },\bar{E})+\overline{h}. & & \displaystyle\end{eqnarray}$$

Inserting $f=f_{\ast }+h$ into the Fokker–Planck equation and using

(8.7) $$\begin{eqnarray}\displaystyle {\dot{f}}_{\ast }=\dot{\unicode[STIX]{x1D713}}_{\ast },\unicode[STIX]{x2202}f_{\ast }/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{\ast }+{\dot{E}}\unicode[STIX]{x2202}f_{\ast }/\unicode[STIX]{x2202}E\simeq c(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701})(\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})+Ze(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t)(\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\bar{E}), & & \displaystyle\end{eqnarray}$$

gives

(8.8) $$\begin{eqnarray}\displaystyle {\dot{h}}-C_{1}\{h-(Mc/Ze)R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\}\simeq (Zef_{M}/T)(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t)-c(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701})(\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}),\quad & & \displaystyle\end{eqnarray}$$

where the lowest-order linearized collision operator must have the property that $C_{1}\{f_{M}\}=0$ . The drive terms for turbulence on the right side are time variation and departure from axisymmetry.

Subtracting out the slowly varying background terms gives the Fokker–Planck equation for the fluctuating response to be

(8.9) $$\begin{eqnarray}\displaystyle \dot{\widetilde{h}} & = & \displaystyle \unicode[STIX]{x2202}\widetilde{h}/\unicode[STIX]{x2202}t+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\widetilde{h}+[\unicode[STIX]{x1D6FA}\boldsymbol{v}\times \boldsymbol{b}-(Ze/M)\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}\tilde{h}\simeq C_{1}\{\tilde{h}\}\nonumber\\ \displaystyle & & \displaystyle +\,(Zef_{M}/T)\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t-c(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701})\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713},\end{eqnarray}$$

with $\unicode[STIX]{x1D6F7}=\bar{\unicode[STIX]{x1D6F7}}+\tilde{\unicode[STIX]{x1D6F7}}$ on the left side and $\tilde{h}=\tilde{h}(\boldsymbol{r},\bar{E},\unicode[STIX]{x1D707},\unicode[STIX]{x1D711},t)$ . Changing to gyrokinetics variables, using $\tilde{h}=\tilde{h}(\boldsymbol{R},\bar{E},\unicode[STIX]{x1D707},t)$ as $\unicode[STIX]{x2202}\tilde{h}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}|_{R}=0$ to lowest order, and then annihilating the next-order equation, yields the desired form of the $\unicode[STIX]{x1D6FF}f=\tilde{f}$ gyrokinetic equation for the non-adiabatic response to be

(8.10) $$\begin{eqnarray}\displaystyle \langle \dot{\widetilde{h}}\rangle _{R}=\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}t}+(v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{M}+\langle \boldsymbol{v}_{E}\rangle _{R})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}=\langle C_{1}\{\tilde{h}\}\rangle _{R}+{\displaystyle \frac{Zef_{M}}{T}}{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}t}}-c{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}},\quad & & \displaystyle\end{eqnarray}$$

where

(8.11) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{E}\rangle _{R}=cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}(\bar{\unicode[STIX]{x1D6F7}}+\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}) & & \displaystyle\end{eqnarray}$$

and

(8.12) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{1}{f_{M}}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}={\displaystyle \frac{1}{p}}{\displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}+{\displaystyle \frac{Ze}{T}}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6F7}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}+\left(\frac{Mv^{2}}{2T}-{\displaystyle \frac{5}{2}}\right){\displaystyle \frac{1}{T}}{\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}. & & \displaystyle\end{eqnarray}$$

In the linear limit, Laplace transforming in time and Fourier decomposing in toroidal angle gives the replacements $\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}/\unicode[STIX]{x2202}t\rightarrow -\text{i}\unicode[STIX]{x1D714}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}$ and $\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}\rightarrow -\text{i}\ell \langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}$ . Defining

(8.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{\ast }^{T}={\displaystyle \frac{c\ell T}{Zef_{M}}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}, & & \displaystyle\end{eqnarray}$$

recovers the usual instability drive $\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\ast }^{T}$ of linear gyrokinetics (Catto Reference Catto1978). Notice that $\unicode[STIX]{x1D714}_{\ast }^{T}$ only depends on the toroidal mode number $\ell$ .

Any neoclassical response $\bar{h}_{p}=\bar{h}_{p}(\unicode[STIX]{x1D713},\bar{E},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})$ in $\bar{h}$ is unimportant as the gyroaverage of

(8.14) $$\begin{eqnarray}\displaystyle \dot{\bar{h}}_{p}=\dot{\bar{E}}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\bar{E}}}+\dot{\unicode[STIX]{x1D707}}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}}+\dot{\unicode[STIX]{x1D713}}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}\simeq -Ze\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\bar{E}}}+\dot{\unicode[STIX]{x1D707}}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}, & & \displaystyle\end{eqnarray}$$

gives

(8.15) $$\begin{eqnarray}\displaystyle \langle \dot{\bar{h}}_{p}\rangle _{R}\simeq -Ze\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}\frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\bar{E}}-ZeB^{-1}\langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}}+\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{R}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{h}_{p}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}\simeq 0, & & \displaystyle\end{eqnarray}$$

since $\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}\simeq \boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\tilde{\unicode[STIX]{x1D6F7}}=-\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}|_{R}$ means $\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}\simeq 0$ , and along with $\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{R}=0$ gives $\langle \dot{\unicode[STIX]{x1D707}}\rangle _{R}=0$ to lowest order.

When a simple local Fourier decomposition of $\tilde{\unicode[STIX]{x1D6F7}}$ is used then

(8.16) $$\begin{eqnarray}\displaystyle \langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}\propto \langle \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}\rangle _{R}=\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}J_{0}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA}), & & \displaystyle\end{eqnarray}$$

giving $\tilde{h}\propto \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}$ . Therefore, when perturbed quasineutrality,

(8.17) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D6F7}}\unicode[STIX]{x1D6F4}{\displaystyle \frac{Z^{2}e^{2}n}{T}}=\unicode[STIX]{x1D6F4}Ze\int \text{d}^{3}v\tilde{h}=\unicode[STIX]{x1D6F4}Ze\int \text{d}^{3}v\langle \tilde{h}\rangle _{r}, & & \displaystyle\end{eqnarray}$$

is evaluated a second set of Bessel functions appears from

(8.18) $$\begin{eqnarray}\displaystyle \langle \tilde{h}\rangle _{r}\propto \langle \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}\rangle _{r}=\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}J_{0}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA}), & & \displaystyle\end{eqnarray}$$

where $\langle \cdots \rangle _{r}$ is the gyroaverage at fixed $\boldsymbol{r}$ .

In a stellarator $\unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D713}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ is used in place of $\unicode[STIX]{x1D713}_{\ast }$ so that

(8.19) $$\begin{eqnarray}\displaystyle f_{\ast }(\unicode[STIX]{x1D6F9},E)\simeq f_{M}(\unicode[STIX]{x1D713},\bar{E})-(Ze\tilde{\unicode[STIX]{x1D6F7}}/T)f_{M}+\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}+\cdots & & \displaystyle\end{eqnarray}$$

and

(8.20) $$\begin{eqnarray}\displaystyle {\dot{f}}_{\ast }=\dot{\unicode[STIX]{x1D6F9}}\unicode[STIX]{x2202}f_{\ast }/\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}+{\dot{E}}\unicode[STIX]{x2202}f_{\ast }/\unicode[STIX]{x2202}E\simeq cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}(\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})+Ze(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}t)(\unicode[STIX]{x2202}f_{M}/\unicode[STIX]{x2202}\bar{E}). & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

As a result, the gyrokinetic equation for a stellarator is

(8.21) $$\begin{eqnarray}\displaystyle \langle \dot{\widetilde{h}}\rangle _{R} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}t}}+(v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{M}+\langle \boldsymbol{v}_{E}\rangle _{R})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}=\langle C_{1}\{\tilde{h}\}\rangle _{R}\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{Zef_{M}}{T}}{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}t}}-\frac{c}{B}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}.\end{eqnarray}$$

In the $\unicode[STIX]{x1D6FF}f$ tokamak and stellarator equations for $\widetilde{h}$ all terms are normally of the same order in a turbulent plasma, and the $cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}$ term in $\langle \boldsymbol{v}_{E}\rangle _{R}$ is the only nonlinearity (Hitchcock & Hazeltine Reference Hitchcock and Hazeltine1978; Frieman & Chen Reference Frieman and Chen1982; Lee Reference Lee1983). This nonlinearity leads to gyroBohm transport as shown in the next section.

9 GyroBohm transport from the fluctuation gyrokinetic equation and zonal flow

In this section the critical balance procedure of Barnes, Parra & Schekochihin (Reference Barnes, Parra and Schekochihin2011) is used to demonstrate that the turbulent transport is gyroBohm within a geometrical factor.

To begin, the eddy size $\unicode[STIX]{x1D6E5}$ is taken as the measure of the typical value of $k_{\bot }^{-1}$ , that is

(9.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\sim k_{\bot }^{-1}, & & \displaystyle\end{eqnarray}$$

and thereby the drives are ordered as

(9.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D714}_{\ast }^{T}\sim \unicode[STIX]{x1D70C}_{i}v_{i}/a\unicode[STIX]{x1D6E5}. & & \displaystyle\end{eqnarray}$$

The adiabatic and non-adiabatic responses are comparable when $\tilde{h}/f_{M}\sim Ze\tilde{\unicode[STIX]{x1D6F7}}/T$ .

According to the mixing length estimate the nonlinear terms become important when $\unicode[STIX]{x1D735}_{\bot }\tilde{h}\sim \unicode[STIX]{x1D735}_{\bot }\,f_{M}$ or

(9.3) $$\begin{eqnarray}\displaystyle \tilde{h}/f_{M}\sim Ze\tilde{\unicode[STIX]{x1D6F7}}/T\sim \unicode[STIX]{x1D6E5}/a. & & \displaystyle\end{eqnarray}$$

Balancing the nonlinear drift term with parallel streaming

(9.4) $$\begin{eqnarray}\displaystyle v_{i}\widetilde{h}/qR\sim v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}\sim (c/B)\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle \widetilde{\unicode[STIX]{x1D6F7}}\rangle _{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}\sim v_{i}\unicode[STIX]{x1D70C}_{i}{\displaystyle \frac{Ze\tilde{\unicode[STIX]{x1D6F7}}\tilde{h}}{T\unicode[STIX]{x1D6E5}^{2}}}, & & \displaystyle\end{eqnarray}$$

then yields the eddy size estimate

(9.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\sim \unicode[STIX]{x1D70C}_{i}(qR/a)\simeq \unicode[STIX]{x1D70C}_{pi}. & & \displaystyle\end{eqnarray}$$

The eddy turnover time $\unicode[STIX]{x1D70F}$ is estimated by balancing the time derivative with the parallel streaming term,

(9.6) $$\begin{eqnarray}\displaystyle \widetilde{h}/\unicode[STIX]{x1D70F}\sim \unicode[STIX]{x2202}\widetilde{h}/\unicode[STIX]{x2202}t\sim v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}\sim v_{i}\widetilde{h}/qR, & & \displaystyle\end{eqnarray}$$

giving

(9.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}\sim qR/v_{i}\sim 1/\unicode[STIX]{x1D714}_{\ast }^{T}. & & \displaystyle\end{eqnarray}$$

The turbulent diffusivity $D_{\text{turb}}\sim \unicode[STIX]{x1D6E5}^{2}/\unicode[STIX]{x1D70F}$ is then estimated to be

(9.8) $$\begin{eqnarray}\displaystyle D_{\text{turb}}\sim (qR/a)\unicode[STIX]{x1D70C}_{i}^{2}v_{i}/a\sim (qR/a)D_{gB}. & & \displaystyle\end{eqnarray}$$

The preceding crude estimate is an aspect ratio larger than gyroBohm. More details with comparisons to ion temperature gradient simulation results are in Barnes et al. (Reference Barnes, Parra and Schekochihin2011). The turbulent levels found by Parra & Barnes (Reference Parra and Barnes2015a ) for ITG momentum diffusivity are also consistent with (9.8).

Nonlinear beating deposits some of the turbulent energy by an inverse cascade into a temporally varying axisymmetric radial mode with a sheared toroidal $\boldsymbol{E}\times \boldsymbol{B}$ drift velocity referred to as a zonal flow. To reach the residual zonal flow level evaluated by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998), the poloidal dependence of the initial condition gives rise to a geodesic acoustic mode (GAM) that damps away. Sugama & Watanabe (Reference Sugama and Watanabe2006) were able to describe this transient behaviour as it damped to the steady state zonal flow level. The residual zonal flow potential has no poloidal variation so no further Landau damping occurs once the transient geodesic acoustic mode (GAM) response to this mode coupling has Landau damped away. The zonal flow shear suppresses turbulence and is responsible for the nonlinear Dimits upshift from the linear threshold for the onset of the ion temperature gradient (ITG) mode. The Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) zonal flow test demonstrated the need for simulations to properly retain ion polarization effects associated with the magnetic drifts. They did this by solving the linear axisymmetric initial value problem associated with the transit average of the gyrokinetic equation

(9.9) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}t}}+(v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{M})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\widetilde{h}={\displaystyle \frac{Zef_{M}}{T}}\frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}t}. & & \displaystyle\end{eqnarray}$$

By specifying poloidally varying gyroaveraged initial conditions they proved that the time asymptotic solution did not damp to zero, but rather relaxed to a residual zonal flow level of

(9.10) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D6F7}}=\tilde{\unicode[STIX]{x1D6F7}}(t=0)/(1+1.64q^{2}/\unicode[STIX]{x1D700}^{1/2}), & & \displaystyle\end{eqnarray}$$

when $k_{\bot }\unicode[STIX]{x1D70C}_{pi}\ll 1$ and $\unicode[STIX]{x1D700}\ll 1$ . Subsequent work considered the role of collisions (Hinton & Rosenbluth Reference Hinton and Rosenbluth1999; Xiao, Catto & Molvig Reference Xiao, Catto and Molvig2007), shaping (Xiao & Catto Reference Xiao and Catto2006), arbitrary $k_{\bot }\unicode[STIX]{x1D70C}$ (Xiao et al. Reference Xiao, Catto and Molvig2007; Monreal et al. Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016), electromagnetic effects (Catto, Parra & Pusztai Reference Catto, Parra and Pusztai2017) and extensions to stellarators (Monreal et al. Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016 and references therein). Xiao’s work is briefly summarized in Xiao, Catto & Dorland (Reference Xiao, Catto and Dorland2007), where comparisons with numerical simulations are presented.

The next section generalizes gyrokinetics to retain electromagnetic effects.

10 Electromagnetic effects in the local gyrokinetic equation

Electromagnetic effects in a local $\unicode[STIX]{x1D6FF}f$ gyrokinetic description are most conveniently handled by writing the fluctuating magnetic field ${\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}$ as

(10.1) $$\begin{eqnarray}\displaystyle {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}=\unicode[STIX]{x1D735}\times ({\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }+\tilde{A}_{\Vert }\boldsymbol{b})\simeq \unicode[STIX]{x1D735}\times {\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }+\unicode[STIX]{x1D735}\tilde{A}_{\Vert }\times \boldsymbol{b}, & & \displaystyle\end{eqnarray}$$

then relating its parallel component to the fluctuating vector potential ${\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}$ by

(10.2) $$\begin{eqnarray}\displaystyle \boldsymbol{b}\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\equiv \tilde{B}_{\Vert }\simeq \boldsymbol{b}\;\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times {\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }. & & \displaystyle\end{eqnarray}$$

In this way, $\tilde{A}_{\Vert }$ and $\tilde{B}_{\Vert }$ turn out to be the two useful dependent variables needed to describe electromagnetic effects. The fluctuating Ampere’s law

(10.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\tilde{A}_{\Vert }\times \boldsymbol{b}+\tilde{B}_{\Vert }\boldsymbol{b})=(4\unicode[STIX]{x03C0}/c)\unicode[STIX]{x1D6F4}Ze\int \text{d}^{3}v\boldsymbol{v}\tilde{f}, & & \displaystyle\end{eqnarray}$$

then gives the parallel component to be

(10.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}\tilde{A}_{\Vert }\simeq -(4\unicode[STIX]{x03C0}/c)\unicode[STIX]{x1D6F4}Ze\int \text{d}^{3}vv_{\Vert }\,\tilde{f}, & & \displaystyle\end{eqnarray}$$

and the perpendicular component as

(10.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\tilde{B}_{\Vert }\times \boldsymbol{b}=(4\unicode[STIX]{x03C0}/c)\unicode[STIX]{x1D6F4}Ze\int \text{d}^{3}v\boldsymbol{v}_{\bot }\,\tilde{f}. & & \displaystyle\end{eqnarray}$$

The gyrokinetic equation now contains

(10.6) $$\begin{eqnarray}\displaystyle [\boldsymbol{E}+c^{-1}\boldsymbol{v}\times (\boldsymbol{B}+{\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\,)]\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}\,f_{\ast }(\unicode[STIX]{x1D713}_{\ast },\bar{E}), & & \displaystyle\end{eqnarray}$$

with

(10.7) $$\begin{eqnarray}\displaystyle \boldsymbol{E}=-\unicode[STIX]{x1D735}(\bar{\unicode[STIX]{x1D6F7}}+\tilde{\unicode[STIX]{x1D6F7}})-c^{-1}\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\,/\unicode[STIX]{x2202}t. & & \displaystyle\end{eqnarray}$$

The $\unicode[STIX]{x1D735}_{v}\bar{E}=M\boldsymbol{v}$ term is easier to treat so is considered first. Removing the adiabatic piece as before leaves the gyroaverage of the fluctuating terms

(10.8) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D712}}}{\unicode[STIX]{x2202}t}}\equiv {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\left(\tilde{\unicode[STIX]{x1D6F7}}-\frac{v_{\Vert }}{c}\tilde{A}_{\Vert }\right)-\frac{1}{c}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\boldsymbol{v}_{\bot }\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}), & & \displaystyle\end{eqnarray}$$

namely,

(10.9) $$\begin{eqnarray}\displaystyle c\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}\equiv c\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}-v_{\Vert }\langle \tilde{A}_{\Vert }\rangle _{R}-\langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\rangle _{R}. & & \displaystyle\end{eqnarray}$$

The gyroaverage $\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}$ can be evaluated when $\tilde{\unicode[STIX]{x1D6F7}}$ , ${\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}$ and ${\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}$ are proportional to $\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}$ . Recalling $\langle \text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}\rangle _{R}=J_{0}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA})$ and using $\langle \boldsymbol{v}_{\bot }\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}\rangle _{R}=\text{i}k_{\bot }^{-1}v_{\bot }\boldsymbol{k}\times \boldsymbol{b}J_{1}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA})$ gives

(10.10) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\rangle _{R}\simeq k_{\bot }^{-1}v_{\bot }J_{1}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA})\tilde{B}_{\Vert }(\boldsymbol{R},t), & & \displaystyle\end{eqnarray}$$

where $\tilde{B}_{\Vert }\simeq \text{i}\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{k}\times {\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }$ . Then using $\tilde{\unicode[STIX]{x1D712}}=\unicode[STIX]{x1D6F4}_{\boldsymbol{k}}\tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}}\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{r}}$ results in the Fourier transform expression

(10.11) $$\begin{eqnarray}\displaystyle \langle \tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}}\rangle _{R}=J_{0}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA})[\tilde{\unicode[STIX]{x1D6F7}}(\boldsymbol{R},t)-c^{-1}v_{\Vert }\tilde{A}_{\Vert }(\boldsymbol{R},t)]+c^{-1}k_{\bot }^{-1}v_{\bot }J_{1}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA})\tilde{B}_{\Vert }(\boldsymbol{R},t),\quad & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D6F7}}(\boldsymbol{R},t)$ , $\tilde{A}_{\Vert }(\boldsymbol{R},t)$ and $\tilde{B}_{\Vert }(\boldsymbol{R},t)$ are proportional to $\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}$ and their $\boldsymbol{k}_{\bot }$ subscripts are suppressed for notational simplicity. To verify that the $k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA}\ll 1$ limit is recovered correctly for the $\tilde{B}_{\Vert }$ term as well as for $\tilde{\unicode[STIX]{x1D6F7}}-c^{-1}v_{\Vert }\tilde{A}_{\Vert }$ , use $\boldsymbol{v}_{\bot }=(\unicode[STIX]{x2202}\boldsymbol{v}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711})\times \boldsymbol{b}$ to integrate by parts, then use $\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }=-\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}|_{R}$ to find

(10.12) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\rangle _{R}=\left\langle \boldsymbol{b}\times \boldsymbol{v}\boldsymbol{\cdot }\left.\frac{\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}\right|_{R}\right\rangle _{R}=-\frac{1}{\unicode[STIX]{x1D6FA}}\langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }\times \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{v}_{\bot }\rangle _{R}\xrightarrow[{k_{\bot }\unicode[STIX]{x1D70C}\ll 1}]{}-\frac{v_{\bot }^{2}\tilde{B}_{\Vert }}{2\unicode[STIX]{x1D6FA}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The preceding results imply the orderings

(10.13) $$\begin{eqnarray}\displaystyle c\tilde{\unicode[STIX]{x1D6F7}}\sim v_{e}\tilde{A}_{\Vert }\sim v_{i}\unicode[STIX]{x1D70C}_{i}\tilde{B}_{\Vert }\sim k_{\bot }\unicode[STIX]{x1D70C}_{i}v_{i}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }, & & \displaystyle\end{eqnarray}$$

giving

(10.14) $$\begin{eqnarray}\displaystyle \tilde{A}_{\Vert }\sim \unicode[STIX]{x1D70C}_{e}\tilde{B}_{\Vert }\sim k_{\bot }\unicode[STIX]{x1D70C}_{e}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }, & & \displaystyle\end{eqnarray}$$

with $v_{e}$ and $v_{i}$ the electron and ion thermal speeds. Recalling $Ze\tilde{\unicode[STIX]{x1D6F7}}/T\sim 1/k_{\bot }a$ , $c\tilde{\unicode[STIX]{x1D6F7}}\sim v_{e}\tilde{A}_{\Vert }$ gives

(10.15) $$\begin{eqnarray}\displaystyle {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}_{\bot }/B\sim k_{\bot }\tilde{A}_{\Vert }/B\sim \unicode[STIX]{x1D70C}_{e}/a, & & \displaystyle\end{eqnarray}$$

and

(10.16) $$\begin{eqnarray}\displaystyle \tilde{B}_{\Vert }/B\sim k_{\bot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }/B\sim 1/k_{\bot }a. & & \displaystyle\end{eqnarray}$$

Consequently, electron dynamics appears to allow electromagnetic effects to enter at very small levels. However, typically only the transit average response of the electrons matters, bringing the electron parallel current down to the ion level. As a result, the parallel Ampere’s law gives

(10.17) $$\begin{eqnarray}\displaystyle k_{\bot }^{2}\tilde{A}_{\Vert }\sim (4\unicode[STIX]{x03C0}\bar{n}_{i}/c)Zev_{i}\tilde{f}/f_{M}, & & \displaystyle\end{eqnarray}$$

which for $k_{\bot }\unicode[STIX]{x1D70C}_{i}$ yields $v_{i}\tilde{A}_{\Vert }/c\sim \unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D6F7}}$ , where $\unicode[STIX]{x1D6FD}$ is the plasma over magnetic pressure, $\unicode[STIX]{x1D6FD}=8\unicode[STIX]{x03C0}nT/B^{2}$ . The transit averaged electron response from $\unicode[STIX]{x2202}\tilde{A}_{\Vert }/\unicode[STIX]{x2202}t$ competes with the electrostatic potential terms in quasineutrality when $\unicode[STIX]{x1D714}\tilde{A}_{\Vert }/k_{\Vert }c\sim \tilde{\unicode[STIX]{x1D6F7}}$ . Combining these two estimates for $\unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D714}_{\ast }\sim k_{\bot }\unicode[STIX]{x1D70C}_{i}v_{i}/a$ and $k_{\Vert }\sim qR$ , suggests that finite $\unicode[STIX]{x1D6FD}$ effects due to shear Alfvén effects enter when

(10.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}=8\unicode[STIX]{x03C0}nT/B^{2}\sim a/qR\ll 1. & & \displaystyle\end{eqnarray}$$

Within other geometrical factors it is very roughly the $\unicode[STIX]{x1D6FD}$ value at which trapped electron modes make the transition to kinetic ballooning modes (Pueschel, Kammerer & Jenko Reference Pueschel, Kammerer and Jenko2008).

To estimate when $\tilde{B}_{\Vert }$ enters, the perpendicular Ampere’s law is used by noting the ion response is what matters in the perpendicular current because of the gyroaverage contained in the integral over velocity space. As a result, it gives

(10.19) $$\begin{eqnarray}\displaystyle k_{\bot }\tilde{B}_{\Vert }\sim (4\unicode[STIX]{x03C0}\bar{n}_{i}/c)Zev_{i}\tilde{f}/f_{M}. & & \displaystyle\end{eqnarray}$$

Then estimating $\tilde{B}_{\Vert }/B\sim 1/k_{\bot }a\sim \tilde{f}/f_{M}$ yields the much higher plasma over magnetic pressure of

(10.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}=8\unicode[STIX]{x03C0}nT/B^{2}\sim 1, & & \displaystyle\end{eqnarray}$$

as the value of $\unicode[STIX]{x1D6FD}$ when compressional Alfvén effects should be kept.

The $\unicode[STIX]{x1D735}_{v}\unicode[STIX]{x1D713}_{\ast }=-(Mc/Ze)R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ term must also be evaluated. However, $\unicode[STIX]{x1D713}\rightarrow \unicode[STIX]{x1D713}+\tilde{\unicode[STIX]{x1D713}}$ with $\tilde{\unicode[STIX]{x1D713}}=-R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}$ must be used in $\unicode[STIX]{x1D713}_{\ast }$ when forming $\dot{\unicode[STIX]{x1D713}}_{\ast }$ so that $\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}/\unicode[STIX]{x2202}t$ cancels the $-c^{-1}\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}/\unicode[STIX]{x2202}t$ in the fluctuating electric field term, leading to the consideration of just

(10.21) $$\begin{eqnarray}\displaystyle R^{2}(\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}-c^{-1}\boldsymbol{v}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}. & & \displaystyle\end{eqnarray}$$

Employing

(10.22) $$\begin{eqnarray}\displaystyle R^{2}v_{\Vert }\boldsymbol{b}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701} & = & \displaystyle -B^{-2}v_{\Vert }\unicode[STIX]{x1D735}\tilde{A}_{\Vert }\boldsymbol{\cdot }\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & = & \displaystyle -B^{-2}v_{\Vert }\unicode[STIX]{x1D735}\tilde{A}_{\Vert }\boldsymbol{\cdot }(I\boldsymbol{B}-R^{2}B^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701})\simeq v_{\Vert }\unicode[STIX]{x2202}\tilde{A}_{\Vert }/\unicode[STIX]{x2202}\unicode[STIX]{x1D701},\end{eqnarray}$$

with the $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{A}_{\Vert }$ term ignored as small, gives

(10.23) $$\begin{eqnarray}\displaystyle R^{2}[\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}-c^{-1}v_{\Vert }\boldsymbol{b}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}={\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}\left(\tilde{\unicode[STIX]{x1D6F7}}-{\displaystyle \frac{v_{\Vert }}{c}}\tilde{A}_{\Vert }\right). & & \displaystyle\end{eqnarray}$$

In addition,

(10.24) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{\bot }\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\rangle _{R}=\langle \boldsymbol{v}_{\bot }\times (\unicode[STIX]{x1D735}\times {\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}})\rangle _{R}=\langle \unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\boldsymbol{\cdot }\boldsymbol{v}_{\bot }\rangle _{R}-\langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }(\unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}})\rangle _{R}\simeq \langle \unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\boldsymbol{\cdot }\boldsymbol{v}_{\bot }\rangle _{R}, & & \displaystyle\end{eqnarray}$$

where the gyroaverage of $\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }=\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }=-\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}_{\bot }/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}|_{R}$ leaves only

(10.25) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{\bot }\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\rangle _{R}\simeq {\displaystyle \frac{1}{R^{2}}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\langle {\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\boldsymbol{\cdot }\boldsymbol{v}_{\bot }\rangle _{R}\simeq -{\displaystyle \frac{v_{\bot }}{R^{2}k_{\bot }}}J_{1}(k_{\bot }v_{\bot }/\unicode[STIX]{x1D6FA}){\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}\tilde{B}_{\Vert }(\boldsymbol{R},t). & & \displaystyle\end{eqnarray}$$

As a result, $\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}$ is again obtained, with

(10.26) $$\begin{eqnarray}\displaystyle \langle R^{2}(\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}-c^{-1}\boldsymbol{v}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\rangle _{R}=\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}. & & \displaystyle\end{eqnarray}$$

Therefore, the right side of the electromagnetic gyrokinetic equation is found to be

(10.27) $$\begin{eqnarray}\displaystyle \langle \dot{\widetilde{h}}\rangle _{R}=C_{1}\{\tilde{h}\}_{R}+{\displaystyle \frac{Zef_{M}}{T}}{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}}{\unicode[STIX]{x2202}t}}-c{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}. & & \displaystyle\end{eqnarray}$$

On the left side of the gyrokinetic equation the electromagnetic nonlinear term must be evaluated, namely

(10.28) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}\left\langle \left(-c\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}-\frac{\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}}{\unicode[STIX]{x2202}t}+\boldsymbol{v}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}\boldsymbol{R}\right\rangle _{R} & = & \displaystyle \boldsymbol{b}\times \left\langle c\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}-\frac{\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D735}({\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\boldsymbol{\cdot }\boldsymbol{v})+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\right\rangle _{R}\nonumber\\ \displaystyle & \simeq & \displaystyle \boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle c\tilde{\unicode[STIX]{x1D6F7}}-v_{\Vert }\tilde{A}_{\Vert }-\boldsymbol{v}_{\bot }\boldsymbol{\cdot }{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\rangle _{R},\end{eqnarray}$$

where the $\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}/\unicode[STIX]{x2202}t$ is small and $\langle \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\rangle _{R}=0$ . As a result, $\langle \boldsymbol{v}_{\bot }\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}\rangle _{R}=\text{i}k_{\bot }^{-1}v_{\bot }\boldsymbol{k}\times \boldsymbol{b}J_{1}$ gives

(10.29) $$\begin{eqnarray}\displaystyle (Ze/M)\langle [-c\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}-\unicode[STIX]{x2202}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}/\unicode[STIX]{x2202}t+\boldsymbol{v}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}\boldsymbol{R}\rangle _{R}\simeq cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}. & & \displaystyle\end{eqnarray}$$

Using the preceding results, the $\unicode[STIX]{x1D6FF}f=\tilde{f}$ form for the nonlinear electromagnetic gyrokinetic equation for the non-adiabatic response is found to be

(10.30) $$\begin{eqnarray}\displaystyle \langle \dot{\widetilde{h}}\rangle _{R} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}t}}+(v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{D}+\langle \boldsymbol{v}_{\unicode[STIX]{x1D712}}\rangle _{R})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}=\langle C_{1}\{\tilde{h}\}\rangle _{R}+{\displaystyle \frac{Zef_{M}}{T}}{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}}{\unicode[STIX]{x2202}t}}\nonumber\\ \displaystyle & & \displaystyle -\,c{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}},\end{eqnarray}$$

with

(10.31) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{D}=\boldsymbol{v}_{M}+cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\bar{\unicode[STIX]{x1D6F7}} & & \displaystyle\end{eqnarray}$$

and

(10.32) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{\unicode[STIX]{x1D712}}\rangle _{R}\equiv cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}=\text{i}cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D6F4}_{\boldsymbol{k}}\boldsymbol{k}\langle \tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}}\rangle _{R}\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{R}}. & & \displaystyle\end{eqnarray}$$

Consequently, to recover the electromagnetic form from the electrostatic form, $\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}$ is simply replaced by $\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}$ .

Employing $\tilde{h}=\unicode[STIX]{x1D6F4}_{\boldsymbol{k}-\boldsymbol{k}^{\prime }}\tilde{h}_{\boldsymbol{k}-\boldsymbol{k}^{\prime }}\text{e}^{\text{i}(\boldsymbol{k}_{\bot }-\boldsymbol{k}_{\bot }^{^{\prime }})\boldsymbol{\cdot }\boldsymbol{R}}$ and $\tilde{\unicode[STIX]{x1D712}}=\unicode[STIX]{x1D6F4}_{\boldsymbol{k}^{\prime }}\tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}^{\prime }}\text{e}^{\text{i}\boldsymbol{k}_{\bot }^{^{\prime }}\boldsymbol{\cdot }\boldsymbol{r}}$ the nonlinear term becomes

(10.33) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{\unicode[STIX]{x1D712}}\rangle _{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}=cB^{-1}\boldsymbol{b}\boldsymbol{\cdot }[\unicode[STIX]{x1D6F4}_{\boldsymbol{k}^{\prime }}\boldsymbol{k}\times \boldsymbol{k}^{\prime }\langle \tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}^{\prime }}\rangle _{R}\tilde{h}_{\boldsymbol{k}-\boldsymbol{k}^{\prime }}]\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{R}}. & & \displaystyle\end{eqnarray}$$

Then the Fourier transformed version of the electromagnetic gyrokinetic equation is

(10.34) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}_{\boldsymbol{k}}}{\unicode[STIX]{x2202}t}}+v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\tilde{h}_{\boldsymbol{k}}+\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{v}_{D}\tilde{h}_{\boldsymbol{k}}+{\displaystyle \frac{c}{B}}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D6F4}_{\boldsymbol{k}^{\prime }}[\boldsymbol{k}_{\bot }\times \boldsymbol{k}_{\bot }^{^{\prime }}\langle \tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}^{\prime }}\rangle _{R}\tilde{h}_{\boldsymbol{k}-\boldsymbol{k}^{\prime }}]\nonumber\\ \displaystyle & & \displaystyle \quad =\langle C_{1}\{\tilde{h}_{\boldsymbol{k}}\}\rangle _{R}+{\displaystyle \frac{Zef_{M}}{T}}\left[\frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{k}}\rangle _{R}}{\unicode[STIX]{x2202}t}+\text{i}\unicode[STIX]{x1D714}_{\ast }^{T}\langle \tilde{\unicode[STIX]{x1D712}}_{\boldsymbol{ k}}\rangle _{R}\right],\end{eqnarray}$$

with $\tilde{h}=\unicode[STIX]{x1D6F4}_{\boldsymbol{k}}\tilde{h}_{\boldsymbol{k}}\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{R}}$ , $\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{r}}=\text{e}^{\text{i}k_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D713}+\text{i}k_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701}-q\unicode[STIX]{x1D717})}$ , $\boldsymbol{k}_{\bot }=k_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}+k_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D701}-q\unicode[STIX]{x1D717})$ , and $\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{k}_{\bot }\times \boldsymbol{k}_{\bot }^{^{\prime }}=B(k_{\unicode[STIX]{x1D713}}^{^{\prime }}k_{\unicode[STIX]{x1D6FC}}-k_{\unicode[STIX]{x1D713}}k_{\unicode[STIX]{x1D6FC}}^{^{\prime }})$ .

Only the $\boldsymbol{k}\times \boldsymbol{b}$ component of the perpendicular Ampere’s law,

(10.35) $$\begin{eqnarray}\displaystyle k_{\bot }^{2}\tilde{B}_{\Vert \boldsymbol{k}}=(4\unicode[STIX]{x03C0}\text{i}/c)\unicode[STIX]{x1D6F4}Ze\boldsymbol{k}\times \boldsymbol{b}\boldsymbol{\cdot }\int {\text{d}^{3}v\langle \boldsymbol{v}_{\bot }\text{e}^{-\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D746}}\rangle }_{r}\tilde{h}_{\boldsymbol{k}}, & & \displaystyle\end{eqnarray}$$

matters, since the adiabatic response does not contribute and both sides of the $\boldsymbol{k}_{\bot }$ component vanish because $\langle \boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{v}_{\bot }\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}\rangle _{r}=0$ and $\unicode[STIX]{x1D735}\tilde{B}_{\Vert }\times \boldsymbol{b}\approx \text{i}\tilde{B}_{\Vert }\boldsymbol{k}\times \boldsymbol{b}$ .

Once again the electromagnetic correction from the neoclassical response $\bar{h}_{p}=\bar{h}_{p}(\unicode[STIX]{x1D713},\bar{E},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E})$ does not matter. The fluctuating ${\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}$ term is proportional to

(10.36) $$\begin{eqnarray}\displaystyle \boldsymbol{v}\times {\displaystyle\mathop{\boldsymbol{B}}\limits_{{\sim}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}\unicode[STIX]{x1D707} & = & \displaystyle B^{-1}v_{\Vert }\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\boldsymbol{b}\times (\unicode[STIX]{x1D735}\times {\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}})\nonumber\\ \displaystyle & \simeq & \displaystyle B^{-1}v_{\Vert }(\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{A}_{\Vert }-\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\boldsymbol{\cdot }\boldsymbol{v}_{\bot })\simeq B^{-1}v_{\Vert }\boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{A}_{\Vert }\end{eqnarray}$$

and is negligible since $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\displaystyle\mathop{\boldsymbol{A}}\limits_{{\sim}}}\boldsymbol{\cdot }\boldsymbol{v}_{\bot }$ is small and $\langle \boldsymbol{v}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{A}_{\Vert }\rangle _{R}=0$ .

Simulations often use $v_{\Vert }$ , rather than $E$ or $\bar{E}$ , as one of the velocity space variables. In such cases, using $Mv_{\Vert }^{2}/2=\bar{E}-Ze\bar{\unicode[STIX]{x1D6F7}}-\unicode[STIX]{x1D707}B$ , the replacement

(10.37) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{R}\widetilde{h}\rightarrow \unicode[STIX]{x1D735}_{R}\widetilde{h}-(Mv_{\Vert })^{-1}(Ze\unicode[STIX]{x1D735}_{R}\bar{\unicode[STIX]{x1D6F7}}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}_{R}B)\unicode[STIX]{x2202}\tilde{h}/\unicode[STIX]{x2202}v_{\Vert }, & & \displaystyle\end{eqnarray}$$

displays the mirror force. In addition, flux tube simulations use the spatial variables $\unicode[STIX]{x1D713}$ ,

(10.38) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D701}-q(\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717}), & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D717}$ or their gyrokinetic counterparts, where now

(10.39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717})={\displaystyle \frac{I(\unicode[STIX]{x1D713})}{q(\unicode[STIX]{x1D713})}}\int _{0}^{\unicode[STIX]{x1D717}}\text{d}\unicode[STIX]{x1D717}^{\prime }/[R^{2}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717}^{\prime })\boldsymbol{B}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717}^{\prime })\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}^{\prime }]. & & \displaystyle\end{eqnarray}$$

The poloidal angle $\unicode[STIX]{x1D717}$ gives the position along the magnetic field. Notice that for this choice

(10.40) $$\begin{eqnarray}\displaystyle \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}=\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}+\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}}=0, & & \displaystyle\end{eqnarray}$$

as desired.

11 Gyrokinetics with flow

The treatment so far retains some flow features by keeping $cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\bar{\unicode[STIX]{x1D6F7}}$ in $\langle \boldsymbol{v}_{E}\rangle _{R}$ . In particular, for sonic flow $\bar{\unicode[STIX]{x1D6F7}}$ is known to be a flux function to lowest order (Hinton & Wong Reference Hinton and Wong1985; Catto, Bernstein & Tessarotto Reference Catto, Bernstein and Tessarotto1987). The order used here is subsonic, but even for diamagnetic flows $\bar{\unicode[STIX]{x1D6F7}}=\text{}\underline{\bar{\unicode[STIX]{x1D6F7}}}(\unicode[STIX]{x1D713})$ to lowest order (Parra et al. Reference Parra, Barnes, Calvo and Catto2012). Consequently, $\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=I\boldsymbol{B}-R^{2}B^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ can be employed to write the small electric drift

(11.1) $$\begin{eqnarray}\displaystyle \langle \text{}\underline{\boldsymbol{v}}_{E}\rangle _{R}\equiv cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}\text{}\underline{\bar{\unicode[STIX]{x1D6F7}}}=(cIB^{-1}\unicode[STIX]{x2202}\text{}\underline{\bar{\unicode[STIX]{x1D6F7}}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})\boldsymbol{b}-(c\unicode[STIX]{x2202}\text{}\underline{\bar{\unicode[STIX]{x1D6F7}}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}, & & \displaystyle\end{eqnarray}$$

as the difference between two $B/B_{p}$ larger flows, one parallel and the other toroidal. As a result, the left side of the gyrokinetic equation can be written as

(11.2) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}\right)+[(v_{\Vert }-IB^{-1}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}})\boldsymbol{b}+\boldsymbol{v}_{M}+\langle {\displaystyle\mathop{\boldsymbol{v}}\limits_{{\sim}}}_{E}\rangle _{R}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}, & & \displaystyle\end{eqnarray}$$

with the toroidal rotation frequency $\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}=\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}(\unicode[STIX]{x1D713})$ defined as

(11.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}=-c\unicode[STIX]{x2202}\text{}\underline{\bar{\unicode[STIX]{x1D6F7}}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713} & & \displaystyle\end{eqnarray}$$

and the fluctuating electric drift

(11.4) $$\begin{eqnarray}\displaystyle \langle {\displaystyle\mathop{\boldsymbol{v}}\limits_{{\sim}}}_{E}\rangle _{R}=cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}_{R}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}. & & \displaystyle\end{eqnarray}$$

The combination

(11.5) $$\begin{eqnarray}\displaystyle D\widetilde{h}/Dt\equiv \unicode[STIX]{x2202}\widetilde{h}/\unicode[STIX]{x2202}t+\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h} & & \displaystyle\end{eqnarray}$$

is the rotating frame time derivative, and redefining $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D701}-q(\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717})-\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}(\unicode[STIX]{x1D713})t$ gives $D\widetilde{h}/Dt=\unicode[STIX]{x2202}\widetilde{h}/\unicode[STIX]{x2202}t|_{\unicode[STIX]{x1D6FC}}$ .

The preceding suggests adding a parallel flow to the Maxwellian to make it a drifting Maxwellian by introducing

(11.6) $$\begin{eqnarray}\displaystyle w_{\Vert }=v_{\Vert }-IB^{-1}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}, & & \displaystyle\end{eqnarray}$$

then in the toroidally rotating frame the curvature drift will include a Coriolis drift as well, since

(11.7) $$\begin{eqnarray}\displaystyle v_{\Vert }^{2}\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b})\simeq (w_{\Vert }^{2}+2w_{\Vert }\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}IB^{-1})\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}), & & \displaystyle\end{eqnarray}$$

with $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}\simeq -\unicode[STIX]{x1D735}\ell nR$ , as in Parra et al. (Reference Parra, Barnes and Catto2011). Continuing to follow their treatment, but using energy rather than parallel velocity as a variable, the energy in the rotating frame is defined as

(11.8) $$\begin{eqnarray}\displaystyle \bar{U}=M[(v_{\Vert }-IB^{-1}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}})^{2}+v_{\bot }^{2}]/2+Ze\bar{\unicode[STIX]{x1D6F7}}-M\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}^{2}R^{2}/2\simeq \bar{E}-Mv_{\Vert }IB^{-1}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}, & & \displaystyle\end{eqnarray}$$

where $M\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}^{2}(R^{2}-I^{2}B^{-2})/2T=M\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}^{2}(R^{2}B_{p}^{2})/2TB^{2}\ll 1$ .

To retain flow with a parallel drifting Maxwellian the replacement $\bar{E}\rightarrow \bar{U}$ is made in $\bar{f}_{\ast }$ to obtain

(11.9) $$\begin{eqnarray}\displaystyle \bar{f}_{\ast }^{\odot }=\unicode[STIX]{x1D702}(\unicode[STIX]{x1D713}_{\ast })[M/2\unicode[STIX]{x03C0}T(\unicode[STIX]{x1D713}_{\ast })]^{3/2}\text{e}^{-\bar{U}/T(\unicode[STIX]{x1D713}_{\ast })} & & \displaystyle\end{eqnarray}$$

and

(11.10) $$\begin{eqnarray}\displaystyle f=f_{\ast }^{\odot }(\unicode[STIX]{x1D713}_{\ast },U)+h, & & \displaystyle\end{eqnarray}$$

with

(11.11) $$\begin{eqnarray}\displaystyle U=M[(v_{\Vert }-IB^{-1}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}})^{2}+v_{\bot }^{2}]/2+Ze\unicode[STIX]{x1D6F7}-M\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}^{2}R^{2}/2\simeq E-Mv_{\Vert }IB^{-1}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}. & & \displaystyle\end{eqnarray}$$

Defining the parallel drifting Maxwellian

(11.12) $$\begin{eqnarray}\displaystyle \bar{f}_{M}^{\odot }(\unicode[STIX]{x1D713},\bar{U}) & = & \displaystyle \bar{f}_{M}^{\odot }=\unicode[STIX]{x1D702}(\unicode[STIX]{x1D713})[M/2\unicode[STIX]{x03C0}T(\unicode[STIX]{x1D713})]^{3/2}\text{e}^{-\bar{U}/T(\unicode[STIX]{x1D713})}\nonumber\\ \displaystyle & \simeq & \displaystyle n(M/2\unicode[STIX]{x03C0}T)^{3/2}\text{e}^{-Mv^{2}/2T+MI\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}v_{\Vert }/TB},\end{eqnarray}$$

where $M\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}^{2}R^{2}B_{p}^{2}/2TB^{2}\ll 1$ is assumed, then

(11.13) $$\begin{eqnarray}\displaystyle \bar{f}_{M}^{\odot }(\unicode[STIX]{x1D713},\bar{U})=f_{M}^{\odot }\simeq \bar{f}_{M}(\unicode[STIX]{x1D713},\bar{E})\text{e}^{MI\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}v_{\Vert }/TB}. & & \displaystyle\end{eqnarray}$$

Observing

(11.14) $$\begin{eqnarray}\displaystyle f_{\ast }^{\odot }(\unicode[STIX]{x1D713}_{\ast },U)\simeq f_{\ast }(\unicode[STIX]{x1D713}_{\ast },E)\text{e}^{MI\unicode[STIX]{x1D6FA}_{\ast }(\unicode[STIX]{x1D713}_{\ast })v_{\Vert }/T(\unicode[STIX]{x1D713}_{\ast })B}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FA}_{\ast }(\unicode[STIX]{x1D713}_{\ast }\rightarrow \unicode[STIX]{x1D713})=\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}$ and $Iv_{\Vert }/B$ a slow function of space, then a new drive term due to flow shear arises when ${\dot{f}}_{\ast }^{\odot }$ is formed since to lowest order

(11.15) $$\begin{eqnarray}\displaystyle {\dot{f}}_{\ast }^{\odot }/f_{\ast }^{\odot }\simeq {\dot{f}}_{\ast }/f_{\ast }+M(Iv_{\Vert }/B)[\unicode[STIX]{x1D6FA}_{\ast }(\unicode[STIX]{x1D713}_{\ast })/T(\unicode[STIX]{x1D713}_{\ast })]^{\cdot }\simeq {\dot{f}}_{\ast }/f_{\ast }+(MIv_{\Vert }/B)\unicode[STIX]{x2202}(\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}/T)/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}.\quad & & \displaystyle\end{eqnarray}$$

The preceding implies that when toroidal flow is retained and the gyrokinetic equation is written in the toroidally rotating frame for $M\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}^{2}R^{2}B_{p}^{2}/2TB^{2}\ll 1$ , the result is

(11.16) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{D\widetilde{h}}{Dt}}+(w_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{M}+\langle {\displaystyle\mathop{\boldsymbol{v}}\limits_{{\sim}}}_{E}\rangle _{R})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h} & = & \displaystyle \langle C_{1}\{\tilde{h}\}\rangle _{R}+{\displaystyle \frac{Zef_{M}^{\odot }}{T}}{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}}{\unicode[STIX]{x2202}t}}\nonumber\\ \displaystyle & & \displaystyle -\,c{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D712}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}}\left(\frac{\unicode[STIX]{x2202}f_{M}^{\odot }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}+{\displaystyle \frac{MIw_{\Vert }f_{M}^{\odot }}{TB}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right),\end{eqnarray}$$

with

(11.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f=\tilde{f}=\tilde{h}-(Ze\tilde{\unicode[STIX]{x1D6F7}}/T)f_{M}^{\odot } & & \displaystyle\end{eqnarray}$$

and

(11.18) $$\begin{eqnarray}\displaystyle \frac{1}{f_{M}^{\odot }}{\displaystyle \frac{\unicode[STIX]{x2202}f_{M}^{\odot }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}={\displaystyle \frac{1}{p}}{\displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}+{\displaystyle \frac{Ze}{T}}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6F7}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}+\left[{\displaystyle \frac{M(w_{\Vert }^{2}+v_{\bot }^{2})}{2T}}-{\displaystyle \frac{5}{2}}\right]{\displaystyle \frac{1}{T}}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$

The preceding form for the gyrokinetic equation requires $1\gg \unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}R/v_{i}\gg B\unicode[STIX]{x1D6FF}/B_{p}\gg \unicode[STIX]{x1D6FF}$ to be valid.

Local gyrokinetic simulations typically solve the gyrokinetic and Maxwell system of equations on a flux tube by specifying $p,T$ and $\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$ and their radial derivatives. Global simulations specify more involved profiles by linking local flux tubes together. As a result, boundary conditions are more complicated for global runs.

Many more details on transforming the gyrokinetic equation to the toroidally rotating frame are given in Parra et al. (Reference Parra, Barnes and Catto2011) where the derivation is performed to higher order to treat intrinsic rotation. Here the purpose is to give some insight into how subsonic toroidal rotation due to a momentum source can be treated when it spins the plasma above the diamagnetic level. In addition, the preceding is a gyrokinetic equation consistent with the important symmetry properties proven and illustrated in Parra, Barnes & Peeters (Reference Parra, Barnes and Peeters2011) and Sugama et al. (Reference Sugama, Watanabe, Nunami and Nishimura2011). These symmetry properties are discussed in appendix A in a more heuristic manner. The consideration of sonic flows results in some simplifications (Able et al. Reference Able, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), but requires a strong source of momentum input. Even when a sonic impurity species is present, a full treatment of momentum transport and profile evolution is required for the non-sonic species to obtain the correct rotation frequency.

12 When are higher-order gyrokinetic variables required?

Local and global simulations are unable to spatially and temporally evolve radial profiles self-consistently. They must be run for specified background profiles (value and derivative) of the plasma density, ion and electron temperatures and radial electric field. If it were not for the lowest flux function that needs to be determined in the electrostatic potential, evolving profiles could be done by linking multiple flux tube simulations using just the conservation of number and energy equations as only particle and heat transport would need to be determined. However, to evaluate the full electrostatic potential requires evaluating the radial flux of toroidal angular momentum in the vicinity of each flux surface and then linking and evolving these surfaces using conservation of toroidal angular momentum. To evaluate the radial flux of toroidal angular momentum and toroidal angular momentum conservation to the requisite order requires a higher-order gyrokinetic description as realized by Catto et al. (Reference Catto, Simakov, Parra and Kagan2008), demonstrated in detail by Parra & Catto (Reference Parra and Catto2008, Reference Parra and Catto2009, Reference Parra and Catto2010a ) and refined further by Parra et al. (Reference Parra, Barnes and Catto2011) and Parra et al. (Reference Parra, Barnes, Calvo and Catto2012). The goal here is the far less ambitious one of demonstrating why higher-order effects are required.

To begin the task, recall that as profiles evolve due to momentum transport, the ambipolarity condition $\langle \boldsymbol{J}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle =0$ must be satisfied on each flux surface, as the displacement current can normally be ignored. Any error resulting in $\langle \boldsymbol{J}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{\text{error}}\neq 0$ acts as a torque on the plasma that results in an unphysical momentum transport satisfying

(12.1) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{c}{V^{\prime }}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}V^{\prime }\unicode[STIX]{x1D6F4}M\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}vf\boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle =\langle \boldsymbol{J}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{\text{error}}\neq 0 & & \displaystyle\end{eqnarray}$$

in the steady state (Parra & Catto Reference Parra and Catto2010b ). Letting $\unicode[STIX]{x1D6F4}M\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}vf\boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle \simeq R^{2}B_{p}\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}$ and $\langle \boldsymbol{J}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{\text{error}}\simeq RB_{p}J_{r}^{\text{error}}$ , implies that the ambipolarity error must be small enough that it and the off-diagonal stress tensor $\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}$ satisfy $\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}/\unicode[STIX]{x2202}r\gg c^{-1}B_{p}J_{r}^{\text{error}}$ or

(12.2) $$\begin{eqnarray}\displaystyle J_{r}^{\text{error}}/env_{i}\ll \unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}\unicode[STIX]{x1D70C}_{pi}/nTa. & & \displaystyle\end{eqnarray}$$

Turbulent momentum transport implies

(12.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}\sim D_{\text{turb}}\unicode[STIX]{x2202}(MnV_{\unicode[STIX]{x1D701}})/\unicode[STIX]{x2202}r, & & \displaystyle\end{eqnarray}$$

where the toroidal flow is diamagnetic with $V_{\unicode[STIX]{x1D701}}\sim v_{i}\unicode[STIX]{x1D70C}_{pi}/a$ . As a result,

(12.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}/nT\sim (\unicode[STIX]{x1D70C}_{pi}/a)^{2}(\unicode[STIX]{x1D70C}_{i}/a), & & \displaystyle\end{eqnarray}$$

and an unphysical torque will be applied by a very small error in the stress tensor or radial current density of

(12.5) $$\begin{eqnarray}\displaystyle J_{r}^{\text{error}}/env_{i}\sim (\unicode[STIX]{x1D70C}_{pi}/a)^{3}(\unicode[STIX]{x1D70C}_{i}/a). & & \displaystyle\end{eqnarray}$$

In addition, a direct evaluation of the off-diagonal stress requires a very high-order evaluation of the fluctuating distribution function

(12.6) $$\begin{eqnarray}\displaystyle \tilde{h}/f_{M}\sim \unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D701}}/nT\sim (\unicode[STIX]{x1D70C}_{pi}/a)^{2}(\unicode[STIX]{x1D70C}_{i}/a). & & \displaystyle\end{eqnarray}$$

Any error in its evaluation will also act as an unphysical torque (Parra & Catto Reference Parra and Catto2010c ). Such an accurate evaluation of $\tilde{h}$ seems unrealistic since the gyrokinetic equation derived herein only gives $\tilde{h}/f_{M}\sim Ze\widetilde{\unicode[STIX]{x1D6F7}}/T\sim 1/k_{\bot }a\sim \unicode[STIX]{x1D70C}_{i}/a$ , while to next-order gyrokinetics gives $\tilde{h}/f_{M}\sim \unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D70C}_{pi}/a^{2}$ . Fortunately, $\tilde{h}/f_{M}\sim \unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D70C}_{pi}/a^{2}$ is sufficient when moments of the exact Fokker–Planck equation are employed to pick up another order in $\unicode[STIX]{x1D70C}_{pi}/a$ , as anticipated by Catto et al. (Reference Catto, Simakov, Parra and Kagan2008) and confirmed more carefully and rigorously by Parra & Catto (Reference Parra and Catto2008, Reference Parra and Catto2010a ). Moments of the full Fokker–Planck equation are routinely employed in neoclassical transport to avoid having to calculate distribution functions to higher order in the gyroradius expansion. And, of course, are the basis of all evaluations of short and long mean free path transport.

The simplest example is radial particle transport as mentioned in the ambipolarity discussion. It is most conveniently evaluated using the coarse grain average of the $R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ component of the electron momentum conservation to find the lowest-order result

(12.7) $$\begin{eqnarray}\displaystyle \langle n_{e}\boldsymbol{V}_{e}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg}=(c/e)\left\langle e{\tilde{n}}_{e}\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}+mR^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}v\boldsymbol{v}C_{ei}\{f_{e}\}\right\rangle _{cg}, & & \displaystyle\end{eqnarray}$$

with $C_{ei}$ the electron–ion collision operator and $m$ the electron mass. For illustrative purposes only electrostatics need be considered. Only the fluctuating electron density ${\tilde{n}}_{e}$ and potential are required to evaluate the turbulent flux. A direct evaluation of $\langle \int \text{d}^{3}vf_{e}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg}$ would require a more accurate $f_{e}$ as the lowest-order neoclassical contribution to $f_{e}$ is odd in $v_{\Vert }$ , the leading classical part of $f_{e}$ also depends on $\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ and thereby results in a vanishing collisional contribution, while $\langle \tilde{f}_{e}\rangle _{cg}=0$ for the turbulent contribution implying that the turbulent correction to $\langle f_{e}\rangle _{cg}\simeq \bar{f}_{e}$ is required. As $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}$ is higher order, the turbulent particle flux is simply

(12.8) $$\begin{eqnarray}\displaystyle c\langle {\tilde{n}}_{e}\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D6F7}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}\rangle _{cg}\simeq \left\langle \int \text{d}^{3}v\tilde{f}_{e}cB^{-1}\boldsymbol{b}\times \unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right\rangle _{cg}. & & \displaystyle\end{eqnarray}$$

Similarly, to evaluate radial heat fluxes the $R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}v^{2}/2$ moment of the full Fokker–Planck equation is also needed. In this case, the lowest-order expression is

(12.9) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{q}_{j}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{cg} & = & \displaystyle \left\langle \int \text{d}^{3}vf_{j}\left(\frac{M_{j}v^{2}}{2}-\frac{5T_{j}}{2}\right)\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right\rangle _{cg}\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{5c}{2Ze}}\left\langle \tilde{p}_{j}\frac{\unicode[STIX]{x2202}\tilde{T}_{j}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\right\rangle _{cg}+\frac{M_{j}}{2}\langle R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\int \text{d}^{3}v\boldsymbol{v}(M_{j}v^{2}-5T_{j})C_{j}\{f_{j}\}\rangle _{cg},\end{eqnarray}$$

where $\tilde{p}_{j}=M_{j}\int \text{d}^{3}vv^{2}\tilde{f}_{j}/3=\bar{n}_{j}\tilde{T}_{j}+{\tilde{n}}_{j}\bar{T}_{j}$ and $\langle \tilde{p}_{j}\unicode[STIX]{x2202}\tilde{T}_{j}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}\rangle _{cg}=(\bar{T}_{j}/\bar{n}_{j})\langle {\tilde{n}}_{j}\unicode[STIX]{x2202}\tilde{p}_{j}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}\rangle _{cg}$ , with $C_{i}=C_{ii}$ for the ions and $C_{e}=C_{ee}+C_{ei}$ for the electrons. Again, this moment of the full kinetic equation obviates the need to evaluate $f_{j}$ to very high order.

The electron particle flux, and electron and ion heat fluxes, require only the lowest-order solutions of the $\unicode[STIX]{x1D6FF}f=\tilde{f}$ gyrokinetic equation for the fluctuations and the lowest-order drift kinetic equation for the collisional modification to the Maxwellian background.

The lowest-order flux function $\bar{\unicode[STIX]{x1D6F7}}\simeq \bar{\unicode[STIX]{x1D6F7}}(\unicode[STIX]{x1D713})$ is required to complete the transport description and determine the complete ion flow and flow shear at a flux surface. It enters particle and heat fluxes but must be evaluated by solving the conservation of toroidal angular momentum equation. The moment procedure becomes more involved for the evaluation of the radial flux of toroidal angular momentum both because a stress tensor is involved and higher-order gyrokinetic variables are required. Even with these complications the moment procedure is essential as it only requires the next-order modifications of gyrokinetics. At present, the next order gyrokinetic treatment combined with a moment approach seems to be the only practical and realistic way to treat momentum transport and profile evolution. Details are beyond the scope of a gyrokinetic tutorial as the mathematics is rather involved since the $R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}\boldsymbol{v}$ and $R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{v}\boldsymbol{v}\boldsymbol{v}$ moments of the full Fokker–Planck equation are required (Parra & Catto Reference Parra and Catto2010b ,Reference Parra and Catto c ), as well as the associated complexity of having to employ higher-order gyrokinetic variables. Nearly all of the required expressions were evaluated by Parra & Catto (Reference Parra and Catto2008, Reference Parra and Catto2009, Reference Parra and Catto2010b ,Reference Parra and Catto c ); Parra et al. (Reference Parra, Barnes and Catto2011); Parra et al. (Reference Parra, Barnes, Calvo and Catto2012); Parra & Barnes (Reference Parra and Barnes2015a ,Reference Parra and Barnes b ); and Calvo & Parra (Reference Calvo and Parra2015); however, the contribution from global features of the radial turbulence profile due to eddy length scale over unperturbed radial scale length corrections remains an active area of research (Parra & Barnes Reference Parra and Barnes2015a ,Reference Parra and Barnes b ) because of the more complicated boundary conditions needed for such a flux tube simulation.

13 Discussion and prospects

This tutorial attempts to summarize gyrokinetics in a manner accessible to interested magnetic fusion practitioners who are not gyrokinetic experts. The orderings appropriate for gyrokinetics and benefits of taking advantage of separation of scale between the slow temporal and spatial background behaviour and the much shorter scales and more rapid time variation of the turbulent fluctuations are highlighted. Then the gyrokinetic equation is derived in its simplest electrostatic forms. The scale separation assumption allows the adiabatic response of each species to be removed from the local gyrokinetic equation for the non-adiabatic response once a lowest-order Maxwellian is shown to be valid because of the existence of pressure or flux surfaces in a tokamak. The electromagnetic and flow modifications to the gyrokinetic equation for the non-adiabatic response are also obtained, along with a brief mention of the stellarator form for the gyrokinetic equation. The discussions stress the separate intrinsic ambipolarity of turbulent and collisional particle transport in an axisymmetric tokamak. Moreover, the development demonstrates that both turbulent and collisional heat fluxes follow from lowest-order gyrokinetic treatments. The inadequacy of lowest-order gyrokinetics only becomes evident when momentum transport and profile evolution need to be treated in a tokamak to determine the unknown flux function in the electrostatic potential. The need to use great care in such cases to avoid introducing errors that torque the plasma is stressed. The prospects of solving a full hybrid fluid (using conservation of number, momentum, and energy equations) and next-order gyrokinetic equation to provide closure is stressed without delving into the challenging and complex formulation that is required. At present this seems to be the most promising path for solving this important problem and using the results to further optimize tokamak performance.

The lowest-order flux function $\bar{\unicode[STIX]{x1D6F7}}\simeq \bar{\unicode[STIX]{x1D6F7}}(\unicode[STIX]{x1D713})$ is required to complete the transport description and determine the complete ion flow and flow shear at a flux surface. It is evaluated by solving the conservation of toroidal angular momentum equation. The moment procedure becomes more involved for the evaluation of the radial flux of toroidal angular momentum both because a stress tensor is involved and higher-order gyrokinetic variables are required. Even with these complications the moment procedure is essential as it only requires the next-order modifications of gyrokinetics. In spite of these subtleties the procedure is more tractable than a full f Hamiltonian method that requires a third-order Hamiltonian – that is, the lowest Hamiltonian plus all corrections up to and including $(\unicode[STIX]{x1D70C}/a)^{3}$ (Parra & Calvo Reference Parra and Calvo2011; Parra et al. Reference Parra, Barnes, Calvo and Catto2012; Calvo & Parra Reference Calvo and Parra2012). At present only the second-order Hamiltonian has been evaluated (Parra & Calvo Reference Parra and Calvo2011; Calvo & Parra Reference Calvo and Parra2015). Therefore, the next-order gyrokinetic treatment combined with a moment approach seems the only practical and realistic way to treat momentum transport and profile evolution.