Skip to main content Accessibility help
×
Home

Information:

  • Access

Actions:

MathJax
MathJax is a JavaScript display engine for mathematics. For more information see http://www.mathjax.org.
      • Send article to Kindle

        To send this article to your Kindle, first ensure no-reply@cambridge.org is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the @free.kindle.com or @kindle.com variations. ‘@free.kindle.com’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘@kindle.com’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Practical gyrokinetics
        Available formats
        ×

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Practical gyrokinetics
        Available formats
        ×

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Practical gyrokinetics
        Available formats
        ×
Export citation

Abstract

Thousands of gyrokinetic papers have been published since the introduction of the gyrokinetic change of variables. The intent here is not to review the field, but rather the goal is to present a tutorial providing insight into why gyrokinetics is an appropriate description of turbulent transport. The focus is on turbulent transport in axisymmetric tokamaks, but many of the ideas and techniques are applicable to stellarators and other magnetic fusion devices with nested surfaces of constant pressure. Besides the origins of gyrokinetics and recent insights, gyrokinetic orderings and gyrokinetic variables are summarized. Then a compact, but careful, derivation of the simplest electrostatic gyrokinetic equation for tokamaks is presented, along with a brief mention of gyrokinetics for stellarators. The advantages of assuming scale separation between the finer spatial scales and faster time variation of the turbulence and the global behaviour and slow temporal evolution of the background are stressed. Moreover, the procedure for removing the adiabatic or Maxwell–Boltzmann response is emphasized. Scale separation allows the near Maxwellian behaviour of the background and the rapid variation of the turbulence to be described by separate, local gyrokinetic equations. The turbulent fluctuations are found by solving the nonlinear gyrokinetic equation for the non-adiabatic portion of the fluctuating distribution function. Also, generalizations of the gyrokinetic equation so that it retains electromagnetic and flow modifications are considered in detail along with some symmetry properties. In addition, ambipolarity, particle transport and heat transport are addressed, along with a discussion of the complications associated with describing momentum transport and profile evolution.

1 Early history

The term gyrokinetics was first used by Catto & Tsang (1977) slightly before the gyrokinetic change of variables was introduced by Catto (1978). Catto & Tsang (1977) 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 (1968) and Taylor & Hastie (1968) that treated small perpendicular wavelengths by using an eikonal technique rather than employing a gyrokinetic change of variables. However, the Catto (1978) 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 (1973) 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 1978), and soon after fully electromagnetic linear descriptions appeared, both without the use of gyrokinetic variables (Antonsen & Lane 1980) and with the use of gyrokinetic variables (Catto, Tang & Baldwin 1981). In addition, nonlinear electromagnetic treatments appeared that employed unperturbed gyrokinetic variables in tokamak geometry (Hitchcock & Hazeltine 1978; Frieman & Chen 1982), and nonlinear electrostatic formulations in slab geometry (Lee 1983) retaining perturbed gyrokinetic variables (Dubin et al. 1983). Hahm, Lee & Brizard (1988) extended the perturbed gyrokinetic variables to nonlinear electromagnetics in a uniform magnetic field $\boldsymbol{B}$ , and Hahm (1988) and Brizard (1989) 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 1982), and the second was the nonlinear gyrokinetic particle-in-cell (PIC) code of Lee (1983) in slab geometry. Years later, Kotschenreuther, Rewoldt & Tang (1995) 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. (2000) 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 1995; Dimits et al. 1996). 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. 2000), 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 (1998) 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 (2006) 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 2003), GENE (Dannert & Jenko 2005), GKV (Sugama & Watanabe 2006) and GKW (Peeters et al. 2009), and even global PIC codes such as ELMFIRE (Heikkinen et al. 2001). Additional global delta f PIC codes like GTC (Lin et al. 2002), GEM (Chen & Parker 2003), GTS (Wang et al. 2006) and ORB5 (Jolliet et al. 2007), 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. (1996) 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 (2009). Appendix A of Sugama & Horton (1998) 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. (2008) and verified by Parra & Catto (2008), as well as the need to satisfy important symmetry properties of the gyrokinetic equation (Parra, Barnes & Catto 2011; Sugama et al. 2011).

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. (2008), Parra & Catto (2008) 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 2010b , 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 2010a , b ; Parra et al. 2011; Parra et al. 2012; Parra & Barnes 2015a , 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 2015a ). Full f descriptions avoid scale separation, but cannot yet be formulated to high enough order in the gyrokinetic variables (Parra et al. 2012). 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 (1992). 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 1973), 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 2008, 2010).

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 (1969) and Rutherford (1970). 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. (1996). 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 (2009). Appendix A of Sugama & Horton (1998) 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 1978). 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 1978; Frieman & Chen 1982; Lee 1983). 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 (2011) 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. (2011). The turbulent levels found by Parra & Barnes (2015a ) 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 (1998), the poloidal dependence of the initial condition gives rise to a geodesic acoustic mode (GAM) that damps away. Sugama & Watanabe (2006) 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 (1998) 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 1999; Xiao, Catto & Molvig 2007), shaping (Xiao & Catto 2006), arbitrary $k_{\bot }\unicode[STIX]{x1D70C}$ (Xiao et al. 2007; Monreal et al. 2016), electromagnetic effects (Catto, Parra & Pusztai 2017) and extensions to stellarators (Monreal et al. 2016 and references therein). Xiao’s work is briefly summarized in Xiao, Catto & Dorland (2007), 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 2008).

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 1985; Catto, Bernstein & Tessarotto 1987). 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. 2012). 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. (2011). 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. (2011) 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 (2011) and Sugama et al. (2011). 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. 2013), 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. (2008), demonstrated in detail by Parra & Catto (2008, 2009, 2010a ) and refined further by Parra et al. (2011) and Parra et al. (2012). 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 2010b ). 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 2010c ). 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. (2008) and confirmed more carefully and rigorously by Parra & Catto (2008, 2010a ). 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 2010b , 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 (2008, 2009, 2010b , c ); Parra et al. (2011); Parra et al. (2012); Parra & Barnes (2015a , b ); and Calvo & Parra (2015); 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 2015a , 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 2011; Parra et al. 2012; Calvo & Parra 2012). At present only the second-order Hamiltonian has been evaluated (Parra & Calvo 2011; Calvo & Parra 2015). 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.

Acknowledgements

P.J.C. is indebted to Professor N. Loureiro and Dr P. Bonoli (both of MIT) and Professor B. Dorland (U. Maryland) for nominating him for a tutorial talk at the 2018 APS Division of Plasma Physics meeting in Portland that served as the impetus for this manuscript; and to Professors F. Parra and M. Barnes (both of Oxford U.) for their support and help with the talk and manuscript. He also appreciates the careful reviewers of JPP whose suggestions and comments further improved the manuscript. The work was supported by the US Department of Energy grant DE-FG02-91ER-54109.

Appendix A. Symmetry properties of the gyrokinetic equation

Parra et al. (2011) and Sugama et al. (2011) have proven that the gyrokinetic equation has important symmetry properties. To simplify the demonstration only an up–down mirror reflected symmetric tokamak is considered for the more heuristic treatment sketched here. For these up–down reflected symmetry considerations the $Z$ axis is the axis of rotation and upon reflection $Z\rightarrow -Z$ , $\unicode[STIX]{x1D735}Z\rightarrow -\unicode[STIX]{x1D735}Z$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\rightarrow -\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ , while $\boldsymbol{B}_{p}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ is unchanged. The low field side of the equatorial plane is taken to be $\unicode[STIX]{x1D717}=0$ . Up-down symmetry means that in the absence of momentum sources or sinks, there can be no radial transport of toroidal angular momentum associated with the lowest-order gyrokinetic equation.

Using $\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]{x1D6FC}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ with $\unicode[STIX]{x1D713}$ , $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D701}-q(\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D717})$ and $\unicode[STIX]{x1D717}$ as the three spatial variables, and considering the linearized electrostatic gyrokinetic equation at first, then the transformations $v_{\Vert }\rightarrow -v_{\Vert }$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\rightarrow -\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ leave the parallel streaming term, as well as the time derivative and collision terms on the left side unchanged. The remaining linear terms are the magnetic drifts

(A 1) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h}=\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}}+\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}}+\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}. & & \displaystyle\end{eqnarray}$$

As the magnetic drift is predominately down (in the $-\unicode[STIX]{x1D735}Z$ direction), up–down symmetry ( $Z\rightarrow -Z$ , $\unicode[STIX]{x1D717}\rightarrow -\unicode[STIX]{x1D717}$ ) means that upon coordinate reflection the magnetic drift becomes predominately in the $\unicode[STIX]{x1D735}Z$ direction with $\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\rightarrow -\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ . As a result, $\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\rightarrow \boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ , but $\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rightarrow -\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ . Moreover, the toroidal drift cannot change direction so $\boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}\rightarrow \boldsymbol{v}_{M}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$ . Therefore, for perturbations the transformations $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}\rightarrow -\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}$ or $k_{\unicode[STIX]{x1D713}}\rightarrow -k_{\unicode[STIX]{x1D713}}$ , $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}\rightarrow \unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}$ , and $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}\rightarrow \unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}$ or $k_{\unicode[STIX]{x1D6FC}}\rightarrow k_{\unicode[STIX]{x1D6FC}}$ do not change the magnetic drift terms and thereby the linearized electrostatic gyrokinetic equation.

When the nonlinear term

(A 2) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{E}\rangle _{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}\widetilde{h} & = & \displaystyle {\displaystyle \frac{c}{B^{2}}}\left\{B^{2}\left[\frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}-{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}}\frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\,\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\left[\frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}}\frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\right]\nonumber\\ \displaystyle & & \displaystyle \left.\quad +\,I\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\left[{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}}\right]\right\}\end{eqnarray}$$

is retained, additional complications arise. As $B^{2}$ is an even function of $\unicode[STIX]{x1D717}$ , a sign change will occur in the first nonlinear term. This observation suggests using the additional transformations $\tilde{h}\rightarrow -\tilde{h}$ and $\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}\rightarrow -\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}$ in all terms of the gyrokinetic equation. Then using $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D717}=R^{-2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}-(\unicode[STIX]{x2202}q/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ , the second nonlinear term transforms properly as well. The third and final nonlinear term appears to be a problem; however, it is an order smaller since it can be rewritten as

(A 3) $$\begin{eqnarray}\displaystyle \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\left[\frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}\frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}-{\displaystyle \frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}}\right]=\frac{\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{h}-\frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}, & & \displaystyle\end{eqnarray}$$

where $cIB^{-1}\unicode[STIX]{x2202}\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}\ll |v_{\Vert }|$ implies that the first term is an order smaller than the parallel streaming term, while in the second term the parallel electric field is small compared to the perpendicular electric field.

The preceding symmetry properties mean that if $\tilde{h}(k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},\bar{E},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E},t)$ and $\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}(k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},t)$ with $\unicode[STIX]{x1D70E}=v_{\Vert }/|v_{\Vert }|$ are solutions to the nonlinear electrostatic gyrokinetic equation then so are $-\tilde{h}(-k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},\bar{E},\unicode[STIX]{x1D707},-\unicode[STIX]{x1D70E},t)$ and $-\langle \tilde{\unicode[STIX]{x1D6F7}}\rangle _{R}(-k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},t)$ . These symmetry properties are then easily extended to the nonlinear electromagnetic gyrokinetic equation by noting that when the solutions $\langle \boldsymbol{A}_{\Vert }\rangle _{R}(k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},t)$ and $\langle \tilde{B}_{\Vert }\rangle _{R}(k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},t)$ are retained, $\langle \boldsymbol{A}_{\Vert }\rangle _{R}(-k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},t)$ and $-\langle \tilde{B}_{\Vert }\rangle _{R}(-k_{\unicode[STIX]{x1D713}},k_{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D717},t)$ are the corresponding solutions. Finally, when subsonic rotation is much stronger than the diamagnetic level, the additional symmetry properties in the rotating frame of $\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}\rightarrow -\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}$ and $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}\rightarrow -\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D701}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$ enter, where $v_{\Vert }$ is replaced by $w_{\Vert }$ .

The direction of the toroidal magnetic field is unimportant for the preceding symmetry argument as the full gyrokinetic equation posses another simpler symmetry property of just $v_{\Vert }\rightarrow -v_{\Vert }$ (or $w_{\Vert }\rightarrow -w_{\Vert }$ ) and $I\rightarrow -I$ .

References

Able, I. G., Plunk, G. G., Wang, E., Barnes, M., Cowley, S. C., Dorland, W. & Schekochihin, A. A. 2013 Multiscale gyrokinetics for rotating tokamak plasmas: fluctuations, transport and energy flows. Rep. Prog. Phys. 76, 116201; 69 pp.
Antonsen, T. M. & Lane, B. 1980 Kinetic equations for low frequency instabilities in inhomogeneous plasmas. Phys. Fluids 23, 12051214.
Barnes, M., Parra, F. I. & Schekochihin, A. A. 2011 Critically balanced ion temperature gradient turbulence in fusion plasmas. Phys. Rev. Lett. 107, 115003115004.
Brizard, A. 1989 Nonlinear gyrokinetic Maxwell-Vlasov equations using magnetic co-ordinates. J. Plasma Phys. 41, 541559.
Calvo, I. & Parra, F. I. 2012 Long-wavelength limit of gyrokinetics in a turbulent tokamak and its intrinsic ambipolarity. Plasma Phys. Control. Fusion 54, 115007, 33 pp.
Calvo, I. & Parra, F. I. 2015 Radial transport of toroidal angular momentum in tokamaks. Plasma Phys. Control. Fusion 57, 075006, 19 pp.
Candy, J. & Waltz, R. E. 2003 An Eulerian gyrokinetic-Maxwell solver. J. Comput. Phys. 186, 545581.
Catto, P. J. 1978 Linearized gyro-kinetics. Plasma Phys. 20, 719722.
Catto, P. J., Bernstein, I. B. & Tessarotto, M. 1987 Ion transport in toroidally rotating tokamak plasmas. Phys. Fluids 30, 27842795.
Catto, P. J., Lee, J.-P. & Ram, A. K. 2017 A quasilinear operator retaining magnetic drift effects in tokamak geometry. J. Plasma Phys. 83, 905830611, 29 pp.
Catto, P. J., Parra, F. I. & Pusztai, I. 2017 Electromagnetic zonal flow residual responses. J. Plasma Phys. 83, 905830402, 38 pp.
Catto, P. J., Simakov, A. N., Parra, F. I. & Kagan, G. 2008 Electrostatic turbulence in tokamaks on transport time scales. Plasma Phys. Control. Fusion 50, 115006, 21 pp.
Catto, P. J., Tang, W. M. & Baldwin, D. E. 1981 Generalized gyrokinetics. Plasma Phys. 23, 639650.
Catto, P. J. & Tsang, K. T. 1977 Linearized gyro-kinetic equation with collisions. Phys. Fluids 20, 396401.
Chen, Y. & Parker, S. E. 2003 A $\unicode[STIX]{x1D6FF}$ f particle method for gyrokinetic simulations with kinetic electrons and electromagnetic perturbations. J. Comput. Phys. 189, 463475.
Dannert, T. & Jenko, F. 2005 Gyrokinetic simulation of collisionless trapped-electron mode turbulence. Phys. Plasmas 12, 072309-8.
Dimits, A. M., Bateman, G., Beer, M. A., Cohen, B. I., Dorland, W., Hammett, G. W., Kim, C., Kinsey, J. E., Kotschenreuther, M., Kritz, A. H. et al. 2000 Comparisons and physics basis of tokamak transport models and turbulence simulations. Phys. Plasmas 7, 969983.
Dimits, A. M., Lodestro, L. L. & Dubin, D. H. E. 1992 Gyroaveraged equations for both the gyrokinetic and drift-kinetic regimes. Phys. Fluids B 4, 271277.
Dimits, A. M., Williams, T. J., Byers, J. A. & Cohen, B. I. 1996 Scalings of ion-temperature-gradient-driven anomalous transport in tokamaks. Phys. Rev. Lett. 77, 7174.
Dorland, W., Jenko, F., Kotschenreuther, M. & Rogers, B. N. 2000 Electron temperature gradient turbulence. Phys. Rev. Lett. 85, 55795582.
Dubin, D. H. E., Krommes, J. A., Oberman, C. & Lee, W. W. 1983 Nonlinear gyrokinetic equations. Phys. Fluids 26, 35243535.
Frieman, E. A. & Chen, L. 1982 Nonlinear gyrokinetic equations for low-frequency electromagnetic waves in general plasma equilibria. Phys. Fluids 25, 502508.
Hahm, T. S., Lee, W. W. & Brizard, A. 1988 Nonlinear gyrokinetic theory for finite-beta plasmas. Phys. Fluids 31, 19401948.
Hahm, T. S. 1988 Nonlinear gyrokinetic equations for tokamak microturbulence. Phys. Fluids 31, 26702673.
Hazeltine, R. D. 1973 Recursive derivation of the drift-kinetic equation. Plasma Phys. 15, 7780.
Heikkinen, J. A., Kiviniemi, T. P., Kurki-Suonio, T., Peeters, A. G. & Sipilä, S. K. 2001 Particle simulation of the neoclassical plasmas. J. Comput. Phys. 173, 527548.
Hinton, F. L. & Rosenbluth, M. N. 1999 Dynamics of axisymmetric (E $\times$ B) and poloidal flows in tokamaks. Plasma Phys. Control. Fusion 41, A653A662.
Hinton, F. L. & Wong, S. K. 1985 Neoclassical ion transport in rotating axisymmetric plasmas. Phys. Fluids 28, 30823098.
Hitchcock, D. A. & Hazeltine, R. D. 1978 Gyrokinetic equation for low and high frequency problems. Plasma Phys. 20, 12411252.
Kagan, G. & Catto, P. J. 2008 Arbitrary poloidal gyroradius effects in tokamak pedestals and transport barriers. Plasma Phys. Control. Fusion 50, 085010, 25 pp.
Kagan, G. & Catto, P. J. 2010 Enhancement of the bootstrap current in a tokamak pedestal. Phys. Rev. Lett. 105, 045002-4.
Kotschenreuther, M., Rewoldt, G. & Tang, W. M. 1995 Comparison of initial value and eigenvalue codes for kinetic toroidal plasma instabilities. Comput. Phys. Commun. 88, 128140.
Kovrizhnykh, L. M. 1969 Transport phenomena in toroidal magnetic systems. Sov. Phys. JETP 29, 475482.
Jolliet, S., Bottino, A., Angelino, P., Hatzky, R., Tran, T. M., McMillan, B. F., Sauter, O., Appert, K., Idomura, Y. & Villard, L. 2007 A global collisionless PIC code in magnetic coordinates. Comp. Phys. Commun. 177, 409425.
Lee, W. W. 1983 Gyrokinetic approach in particle simulation. Phys. Fluids 26, 556562.
Lee, X. S., Myra, J. R. & Catto, P. J. 1983 General frequency gyrokinetics. Phys. Fluids 26, 223229.
Lin, Z., Ethier, S., Hahm, T. S. & Tang, W. M. 2002 Size scaling of turbulent transport in magnetically confined plasmas. Phys. Rev. Lett. 88, 195004-4.
Monreal, P., Calvo, I., Sánchez, E., Parra, F. I., Bustos, A., Könies, A., Kleiber, R. & Görler, T. 2016 Residual zonal flows in tokamaks and stellarators at arbitrary wavelengths. Plasma Phys. Control. Fusion 58, 045018, 15 pp.
Parra, F. I. & Barnes, M. 2015a Intrinsic rotation in tokamaks: theory. Plasma Phys. Control. Fusion 57, 045002, 45 pp.
Parra, F. I. & Barnes, M. 2015b Equivalence of two different approaches to global $\unicode[STIX]{x1D6FF}$ f gyrokinetic simulations. Plasma Phys. Control. Fusion 57, 054003, 11 pp.
Parra, F. I., Barnes, M., Calvo, I. & Catto, P. J. 2012 Intrinsic rotation with gyrokinetic models. Phys. Plasmas 19, 056116, 23 pp.
Parra, F. I., Barnes, M. & Catto, P. J. 2011 Sources of intrinsic rotation in the low-flow ordering. Nucl. Fusion 51, 113001, 10 pp.
Parra, F. I., Barnes, M. & Peeters, A. G. 2011 Up-down symmetry of the turbulent transport of toroidal angular momentum in tokamaks. Phys. Plasmas 18, 062501-14.
Parra, F. I. & Calvo, I. 2011 Phase-space Lagrangian derivation of electrostatic gyrokinetics in general geometry. Plasma Phys. Control. Fusion 53, 045001, 49 pp.
Parra, F. I. & Catto, P. J. 2008 Limitations of gyrokinetics on transport time scales. Plasma Phys. Control. Fusion 50, 065014, 23 pp.
Parra, F. I. & Catto, P. J. 2009 Vorticity and intrinsic ambipolarity in turbulent tokamaks. Plasma Phys. Control. Fusion 51, 095008, 38 pp.
Parra, F. I. & Catto, P. J. 2010a Turbulent transport of toroidal angular momentum in low flow gyrokinetics. Plasma Phys. Control. Fusion 52, 045004, 32 pp; Erratum 2010 Plasma Phys. Control. Fusion 52, 059801.
Parra, F. I. & Catto, P. J. 2010b Transport of momentum in full f gyrokinetics. Phys. Plasmas 17, 056106-11.
Parra, F. I. & Catto, P. J. 2010c Non-physical momentum sources in slab geometry gyrokinetics. Plasma Phys. Control. Fusion 52, 085011, 12 pp.
Peeters, A. G., Strintzi, D., Camenen, Y., Angioni, C., Casson, F. J., Hornsby, W. A. & Snodin, A. P. 2009 Influence of the centrifugal force and parallel dynamics on the toroidal momentum transport due to small scale turbulence in a tokamak. Phys. Plasmas 16, 042310-10.
Pueschel, M. J., Kammerer, M. & Jenko, F. 2008 Gyrokinetic turbulence simulations at high plasma beta. Phys. Plasmas 15, 102310, 10 pp.
Rewoldt, G., Tang, W. M. & Chance, M. S. 1982 Electromagnetic kinetic toroidal eigenmodes for general magnetohydrodynamicequilibria. Phys. Fluids 25, 480490.
Rosenbluth, M. N. & Hinton, F. L. 1998 Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks. Phys. Rev. Lett. 80, 724727.
Rutherford, P. H. 1970 Collisional diffusion in an axisymmetric torus. Phys. Fluids 13, 482489.
Rutherford, P. H. & Frieman, E. A. 1968 Drift instabilities in general magnetic field configurations. Phys. Fluids 11, 569585.
Sugama, H., Okamoto, M., Horton, W. & Wakatani, M. 1996 Transport processes and entropy production in toroidal plasmas with gyrokinetic electromagnetic turbulence. Phys. Plasmas 3, 23792394.
Sugama, H. & Horton, W. 1998 Nonlinear electromagnetic gyrokinetic equation for plasma with large mean flows. Plasma Physics 5, 25602573.
Sugama, H. & Watanabe, T.-H. 2006 Collisionless damping of geodesic acoustic modes. J. Plasma Phys. 72, 825828.
Sugama, H., Watanabe, T. H., Nunami, M. & Nishimura, S. 2011 Momentum balance and radial electric fields in axisymmetric and nonaxisymmetric toroidal plasmas. Plasma Phys. Control. Fusion 53, 024004, 17 pp.
Sydora, R. D. 1995 Toroidal gyrokinetic particle simulations of core fluctuations and transport. Phys. Scr. 52, 474480.
Taylor, J. B. & Hastie, R. J. 1968 Stability of general plasma equilibria - I Formal theory. Plasma Phys. 10, 479494.
Wang, W. X., Lin, Z., Tang, W. M., Lee, W. W., Ethier, S., Lewandowski, J. L. V., Rewoldt, G., Hahm, T. S. & Manickam, J. 2006 Gyro-kinetic simulation of global turbulent transport properties in tokamak experiments. Phys. Plasmas 13, 092505-12.
Xiao, Y. & Catto, P. J. 2006 Plasma shaping effects on the collisionless residual zonal flow level. Phys. Plasmas 13, 082307-5.
Xiao, Y. & Catto, P. J. 2006 Short wavelength effects on the collisionless neoclassical polarization and residual zonal flow level. Phys. Plasmas 13, 102311-11.
Xiao, Y., Catto, P. J. & Dorland, W. 2007 Effects of finite poloidal gyroradius, shaping, and collisions on the zonal flow residuals. Phys. Plasmas 14, 055910-6.
Xiao, Y., Catto, P. J. & Molvig, K. 2007 Collisional damping for ion temperature gradient mode driven zonal flow. Phys. Plasmas 14, 032302-11.