2 Background

The quasilinear operator must be valid for a wide range of wave frequencies $\unicode[STIX]{x1D714}$ . For waves harmonic in time and space, our orderings must allow us to consider low frequencies of order

to recover the unmagnetized and low frequency limit of Kennel & Engelmann (Reference Kennel and Engelmann1966). Here and hereafter, $k_{\Vert }$ and $\boldsymbol{k}_{\bot }$ are the parallel and perpendicular components of the wave vector, and $v_{\Vert }$ and $\boldsymbol{v}_{d}$ are the parallel and magnetic drift velocities. For higher frequencies we must allow

to recover magnetized Kennel & Engelmann (Reference Kennel and Engelmann1966), where $\unicode[STIX]{x1D6FA}$ is the gyrofrequency. The large range of frequencies allows our quasilinear result to treat frequencies from drift waves to ion and then electron cyclotron frequencies and above.

For magnetic drifts to modify the usual quasilinear (QL) resonance condition $\unicode[STIX]{x1D714}-p\unicode[STIX]{x1D6FA}-k_{\Vert }v_{\Vert }=0$ requires $\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{v}_{d}\sim k_{\Vert }v_{\Vert }$ , where $p$ is an integer or zero. Letting $\unicode[STIX]{x1D70C}$ and $R$ denote the gyroradius and major radius of the tokamak, we find that

with $\unicode[STIX]{x1D70C}_{p}=\unicode[STIX]{x1D70C}B/B_{p}$ the poloidal gyroradius and $B_{p}$ the poloidal magnetic field. Typically,

so we expect magnetic drift effects will be of most interest when

For large values of $k_{\bot }\unicode[STIX]{x1D70C}_{i}$ we expect the QL diffusivity to be small since it depends on Bessel functions that for large argument behave as $\text{J}_{p}(w)\sim w^{-1/2}$ , with $w\sim k_{\bot }\unicode[STIX]{x1D70C}_{i}$ . As a result, the magnetic drift corrections we evaluate may be important when $k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$ for $B_{p}\ll B$ . Indeed, this seems to explain why the Kennel & Engelmann (Reference Kennel and Engelmann1966) limit of the QL operator works so well in a tokamak.

The quasilinear diffusivity in the constant magnetic field case is proportional to a delta function whose argument is $\unicode[STIX]{x1D714}-p\unicode[STIX]{x1D6FA}-k_{\Vert }v_{\Vert }$ . We desire to generalize this QL diffusivity to include magnetic drifts. We expect to find that generalized diffusivity contains a drift modified resonance condition so the argument of the delta function depends on $\unicode[STIX]{x1D714}-p\unicode[STIX]{x1D6FA}-k_{\Vert }v_{\Vert }-\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{v}_{d}$ . This new resonance condition is now what must be used to gather terms in the proper way to obtain the required form of the quasilinear diffusivity with magnetic drift effects retained.

3 Formulation of the quasilinear and linear equations

Consider the full kinetic equation

with $\boldsymbol{e}$ and $\boldsymbol{b}$ the applied radio frequency (rf) wave fields and $\boldsymbol{B}=B\boldsymbol{n}$ the unperturbed magnetic field. We then take

with

and $f_{0}$ gyrophase independent so that

with $\unicode[STIX]{x1D711}$ the gyrophase. We are interested in situations for which collisions and quasilinear diffusion are equally important in the evolution of $f_{0}$ . We also require that $f_{0}$ does not have any fast time or space dependence so take it to be the coarse grain average of $f$ :

where the integral over toroidal angle $\unicode[STIX]{x1D701}$ insures axisymmetry, and the average $\langle \ldots \rangle _{\unicode[STIX]{x1D717}}$ over poloidal angle $\unicode[STIX]{x1D717}$ removes the high poloidal mode number variation. The integral over $\unicode[STIX]{x1D70F}$ removes fast time variation by taking $\unicode[STIX]{x1D714}^{-1}\ll T\ll \unicode[STIX]{x1D708}^{-1}$ , and the integral over the poloidal flux function $\unicode[STIX]{x1D713}$ removes rapid radial variation by assuming $\unicode[STIX]{x1D706}_{rf}^{\text{radial}}\ll \unicode[STIX]{x1D6E5}\ll a$ , with $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D714}$ being the collision frequency and typical rf wave frequency, respectively, and a and $\unicode[STIX]{x1D706}_{rf}^{\text{radial}}$ being the minor radius and the typical radial wavelength of the rf wave, respectively. More details on how $\langle \ldots \rangle _{\unicode[STIX]{x1D717}}$ deals with rapid variation will be given in § 5, but for now it is sufficient to know that $f_{0}$ is allowed to be a slow function of $\unicode[STIX]{x1D717}$ , as well as $\unicode[STIX]{x1D713}$ and $t$ . Consequently, we are assuming

with $\unicode[STIX]{x1D6FA}=ZeB/Mc$ and $Q$ the quasilinear operator to be derived. We allow $f_{1}$ to rapidly vary in time and space by permitting

These assumptions mean that $f_{0}$ satisfies the quasilinear (QL) equation

while $f_{1}$ satisfies the linear equation

with collisions neglected since we are normally interested in wave frequencies much higher than the collision frequencies. If the sum of these last two equations is subtracted from the full kinetic equation then the kinetic equation for the difference $f-(f_{0}+f_{1})$ is recovered. However, this difference equation is not needed to derive the quasilinear diffusivity. Except for tokamak geometry, equations (3.8) and (3.9) are the standard starting equations for all quasilinear treatments of rf heating and current drive.

For the QL equation we use the drift kinetic limit of gyrokinetics in which we keep the magnetic moment to higher order. In these $\boldsymbol{r},v,\unicode[STIX]{x1D707},\unicode[STIX]{x1D711},t$ drift kinetic variables the magnetic moment is given by

with

the usual lowest-order magnetic moment. We write the velocity $\boldsymbol{v}$ in the usual cylindrical velocity space variables aligned with the magnetic field

with $\unicode[STIX]{x1D711}$ the gyrophase,

the magnitude of velocity of the charge (a constant of the motion for our treatment), and the three orthonormal unit vectors satisfying

Then the next-order correction to the magnetic moment is (Catto et al. Reference Catto, Tang and Baldwin1981; Lee et al. Reference Lee, Myra and Catto1983; Kagan & Catto Reference Kagan and Catto2008; Parra & Catto Reference Parra and Catto2008)

with

the magnetic drift velocity. The unperturbed $\boldsymbol{E}\times \boldsymbol{B}$ drift can be added to $\boldsymbol{v}_{d}$ if the unperturbed electric field $\boldsymbol{E}$ is known. In typical quasilinear treatments, including the treatment here, the perturbed drift due to the applied wave fields is considered to be small and neglected. However, these and other small fluctuating corrections to the unperturbed gyrokinetic trajectories are expected to provide the stochasticity necessary to randomize kicks as charges make successive passes through a resonance (Becoulet, Gambier & Samain Reference Becoulet, Gambier and Samain1991). While low frequency gyrokinetic treatments have been very successful in retaining the effects of low frequency fluctuating fields (Kagan & Catto Reference Kagan and Catto2008; Parra & Catto Reference Parra and Catto2008; Parra & Calvo Reference Parra and Calvo2011; Calvo & Parra Reference Calvo and Parra2012), at present there seems to be no systematic way of treating high frequency fluctuations in high frequency gyrokinetics.

Changing to these new variables, gyroaveraging, and coarse grain averaging results in the QL equation in its drift kinetic form

with the QL operator given by

and where the perturbed acceleration is defined as

There are two important points to mention here. The first point is that $Q\{\,f_{0}\}$ must be allowed to be a slow function of $\unicode[STIX]{x1D717}$ in a tokamak so that $\langle \ldots \rangle _{\unicode[STIX]{x1D717}}\neq (2\unicode[STIX]{x03C0})^{-1}\oint \text{d}\unicode[STIX]{x1D717}(\ldots )$ . The second point is that the drift kinetic treatment requires that we average the QL term over gyrophase, as well as perform a coarse grain average. Therefore, in the unperturbed distribution function

$\unicode[STIX]{x1D70E}=v_{\Vert }/|v_{\Vert }|$ is the sign of $v_{\Vert }$ , $\unicode[STIX]{x1D713}$ is the poloidal flux function and $\unicode[STIX]{x1D717}$ is the poloidal angle variable. The preceding form and/or the original full $f$ form of the kinetic equation are convenient for forming conservation equations, as well as other moment equations. Which form is to be used depends on the accuracy required.

It is important to mention here that we will not obtain a radial spatial derivative in our quasilinear operator since we will evaluate $Q\{\,f_{0}\}$ in drift kinetic variables by assuming

with $\unicode[STIX]{x1D70C}_{p}=\unicode[STIX]{x1D70C}B/B_{p}$ the poloidal gyroradius and $B_{p}$ the poloidal magnetic field.

To evaluate the lowest-order distribution function it is sometimes convenient to solve the transit averaged version of the preceding equation. To lowest order we then assume that streaming dominates so that

with the lowest-order distribution function defined to be independent of poloidal angle

Then

Letting

with

the next-order QL equation becomes

It is then convenient to introduce the transit average

with $\text{d}\unicode[STIX]{x1D70F}=\text{d}\unicode[STIX]{x1D717}B/v_{\Vert }\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}$ and $\bar{f}_{0}$ the transit average of $f_{0}$ . The $\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D717}$ integrations in both the numerators and denominators are over a full poloidal circuit following a charged particle. In this way $v_{\Vert }$ and $\unicode[STIX]{x1D717}$ change signs together at a turning point for trapped particles, and odd functions of $v_{\Vert }$ result in a vanishing transit average (for example, $\overline{v_{\Vert }}=0$ and $\overline{v_{\Vert }/B}=0$ ).

Equation (3.27) treats magnetic drift departures from a flux surface as being small by assuming $\unicode[STIX]{x1D6FF}f_{0}\sim (Iv_{\Vert }/\unicode[STIX]{x1D6FA})\unicode[STIX]{x2202}\bar{f}_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$ so that poloidal gyroradius corrections are small compared to the radial scale length associated with $\bar{f}_{0}$ . This small poloidal gyroradius assumption means that finite orbit effects will be treated as being negligible.

Transit averaging the equation for $\unicode[STIX]{x1D6FF}f_{0}$ and using the periodicity of $\bar{f}_{0}$ we obtain

where $\unicode[STIX]{x2202}\bar{f}_{0}/\unicode[STIX]{x2202}t$ is often ignored as small. This transit averaged equation is the QL equation that must be solved numerically for transit averaged $\bar{f}_{0}$ once the QL term is made explicit.

To solve the linear equation we use the gyrokinetic variables $\boldsymbol{R},v,\unicode[STIX]{x1D707},\unicode[STIX]{x1D711}$ and $t$ of Lee et al. (Reference Lee, Myra and Catto1983) and Parra & Catto (Reference Parra and Catto2008), where

and $v=|\boldsymbol{v}|,\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D707}_{1}$ , $\boldsymbol{v}=v_{\Vert }\boldsymbol{n}+v_{\bot }(\boldsymbol{e}_{1}\cos \unicode[STIX]{x1D711}+\boldsymbol{e}_{2}\sin \unicode[STIX]{x1D711})$ and $\boldsymbol{e}_{1}\times \boldsymbol{e}_{2}=\boldsymbol{n}=\boldsymbol{B}/B$ defined by (3.10) to (3.14), that is, exactly as for drift kinetics.

Higher-order terms not shown in the gyrokinetic variable $\boldsymbol{R}$ are needed to remove the gyrophase dependence from $\dot{\boldsymbol{R}}$ to the requisite order and thereby obtain

where for an arbitrary scalar or vector function $X$ ,

Similarly, the $\unicode[STIX]{x1D707}_{1}$ modification of $\unicode[STIX]{x1D707}_{0}$ results in

to the requisite order. In addition, we can even retain the gyrophase to higher order,

where to lowest order

while to next order

We assume an axisymmetric unperturbed magnetic field of the form

with $\unicode[STIX]{x1D701}$ the toroidal angle variable and $I=RB_{t}$ a flux function, where $B_{t}$ is the toroidal magnetic field. We are free to let $\boldsymbol{e}_{1}$ be the unit vector in the $\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ direction by taking

with $B_{p}$ the poloidal magnetic field and $\boldsymbol{e}_{2}=\boldsymbol{n}\times \boldsymbol{e}_{1}$ .

Using these higher-order gyrokinetic variables the perturbed kinetic equation becomes the high frequency gyrokinetic equation (Lee et al. Reference Lee, Myra and Catto1983)

to the requisite order. To obtain the preceding it is not necessary to gyroaverage the perturbed kinetic equation so that $f_{1}$ depends on guiding centre location $\boldsymbol{R}$ ,

but the rf electric (and magnetic) field depends only on particle location $\boldsymbol{r}$ ,

This distinction between the guiding centre location $\boldsymbol{R}$ and the particle location $\boldsymbol{r}$ gives rise to magnetic drift and Bessel function modifications in (3.39).

4 Gyrokinetic solution of the linear equation

There is an important subtlety we must address before proceeding any further. As we need to perform the gyrokinetic change of variables on the perturbed equation to higher order than usual it is necessary to depart from the more straightforward treatment by solving explicitly for certain pieces of the perturbed distribution function. To do so we first use Faraday’s law,

to write

and then use

to write

We begin by ignoring the poloidal gyroradius over unperturbed scale length corrections by using the approximation

The preceding term requires special handling since it must be treated in the same way we treat the adiabatic portion of the $\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{1}$ term in standard gyrokinetics.

To deal with time derivatives without Fourier (or Laplace) transforming in time we form the time derivative of the perturbed equation

The next step is unconventional, but is one of two steps that allow us to obtain results with drifts included that would not otherwise be possible. We extract the first non-resonant portion of the perturbed distribution function by defining $h$ via

to obtain the alternate and convenient form

where we define

In the preceding equation for $h$ we use

where

with

Next, we use the second unconventional step to extract a second term by using

where we neglect $k_{\bot }^{-1}\unicode[STIX]{x1D735}\ell nB\sim 1/k_{\bot }R\ll 1$ corrections. Then, we extract the second non-resonant piece of the perturbed distribution function by defining $g$ as

to obtain the more convenient form

with

Three of the terms depending on $\boldsymbol{G}\sim \boldsymbol{v}_{d}$ can be seen to be small since

and

Consequently, we are left to consider

We can Fourier decompose the applied electric field $\boldsymbol{e}$ poloidally and toroidally. We also Fourier decompose in time to permit more than a single wave frequency as might be the case for lower hybrid and ion cyclotron heating. In addition, we allow the local fine scale radial eikonal variation to have multiple roots (that must be summed over) due to mode conversion. Therefore, we take

To solve for $f_{1}$ we seek solutions of the form

and

with

and

and where we define

and

The sums are over all wavenumbers and frequencies for the poloidally periodic applied wave fields. However, because of the presence of trapped particles, $g,h$ and $f_{1}$ are not Fourier decomposed poloidally. The extra sum over $p$ in $g,h$ and $f_{1}$ is for the Fourier decomposition in gyrophase. The toroidal and poloidal wave numbers are taken to be $n$ and $m$ , respectively, and $\unicode[STIX]{x1D714}$ denotes wave frequencies. The vector

is the radial wave vector.

The distinction between drift kinetic and gyrokinetic variables is important and must be treated carefully in exponentials. In particular, we need to use the preceding expressions to write the right-hand side of the $g$ equation in terms of the gyrokinetic variables $\unicode[STIX]{x1D6F9},\unicode[STIX]{x1D6E9},\text{Z},v,\unicode[STIX]{x1D707},\unicode[STIX]{x1D711}$ . Before doing so it is useful to note that $f_{p}$ and $g_{p}$ are related by

Using the preceding expressions, multiplying by $\exp (\text{i}p^{\prime }\unicode[STIX]{x1D711})$ , integrating over $\unicode[STIX]{x1D719}$ from 0 to 2 $\unicode[STIX]{x03C0}$ , defining,

and

recalling the ‘Bessel generating function’

and performing the gyrophase integrals using

we obtain in gyrokinetic variables that

where we define

and we have ignored small poloidal gyroradius corrections that arise from Taylor expanding $B^{-1}v_{\Vert }\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D707}$ about $\unicode[STIX]{x1D6F9}$ . The distinction between drift and gyrokinetic variables on the right- and left-hand sides of the equation results in the term $L=\unicode[STIX]{x1D6FA}^{-1}\boldsymbol{v}\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{k}\sim k_{\bot }\unicode[STIX]{x1D70C}$ .

We are now ready to Fourier transform the kinetic equation for $g$ in $\unicode[STIX]{x1D711}$ to obtain

where we define

and

In writing down this form of the kinetic equation for $g_{p}$ we have taken the liberty of using the eikonal approximation for the radial variation to streamline the solution. We could avoid doing this by performing a more complicated trajectory integral in $\unicode[STIX]{x1D6F9}$ as well as $\unicode[STIX]{x1D6E9}$ , rather than just in $\unicode[STIX]{x1D6E9}$ .

We again multiply by $\exp (\text{i}p^{\prime }\unicode[STIX]{x1D711})$ and integrate over $\unicode[STIX]{x1D711}$ , but this time we use the additional results

and

where we define $\boldsymbol{e}_{\pm }$ as

and $\boldsymbol{e}_{\bot p}$ as

with

The alternate form of $\boldsymbol{e}_{\bot p}$ is found by using the Bessel recurrence relations.

Using the preceding we obtain

The approximate sign in (4.50) is a reminder of the approximations (4.5) and (4.13) used to remove the non-resonant terms. An alternate form follows using $q\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}=\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}=I/R^{2}$ with $q=q(\unicode[STIX]{x1D713})$ the safety factor,

where

and

To see that it does not vanish, the integrals needed to evaluate $\unicode[STIX]{x1D6FF}_{p,m}$ can be performed as sketched in the appendix. However, there is actually no need to do so since the $\unicode[STIX]{x1D6FF}_{p,m}$ term can be neglected for the following two reasons. First, to derive the quasilinear operator we will show that only the residual associated with $\unicode[STIX]{x1D714}-k_{\Vert }v_{\Vert }-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}=p\unicode[STIX]{x1D6FA}$ is needed for $p\neq 0$ . Second, for p $=$ 0 we need only assume $\unicode[STIX]{x1D6E5}_{m}\ll \unicode[STIX]{x1D714}.$ The first condition means we can neglect the $\unicode[STIX]{x1D6FF}_{p,m}$ for $p\neq 0$ as long as we assume

The second condition for neglecting $\unicode[STIX]{x1D6FF}_{p,m}$ requires the even less stringent assumption

for $p=0$ , where $\unicode[STIX]{x1D714}_{\ast }\sim k_{\bot }\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x1D6FA}/L_{n}$ is the diamagnetic drift frequency with $L_{n}$ the density scale length.

In addition, we are continuing to assume that on the right-hand side of the kinetic equation for $g_{p}$ the distinction between gyrokinetic and drift kinetic variables is unimportant for $W_{p,m}$ as long as poloidal gyroradius over unperturbed radial scale lengths are small as is normally assumed in the core. Consequently, on the right we may assume $\boldsymbol{R}\simeq \boldsymbol{r}$ to the order required except in $\exp (\text{i}m\unicode[STIX]{x1D6E9})$ . This $\unicode[STIX]{x1D70C}_{p}f_{0}^{-1}\unicode[STIX]{x1D735}f_{0}\ll 1$ assumption may fail in a high confinement pedestal or for the very energetic ions produced by minority heating, but such considerations are beyond the scope of the present treatment.

To ignore the $\unicode[STIX]{x1D70C}/R$ terms in $\bar{\unicode[STIX]{x1D6FA}}$ on the left-hand side of (4.51) for $g_{p}$ , we need only recall that the poloidal variation of the magnetic field is much more important because $B\propto R^{-1}\propto [1+(r/R)\cos \unicode[STIX]{x1D717}]^{-1}$ . As a result, we need only assume

with $\unicode[STIX]{x1D70C}_{i}$ the ion gyroradius. As this assumption is easily satisfied from here on we use

The preceding assumptions mean we need only solve

To integrate the $g_{p}$ equation we introduce the trajectory time variable $\unicode[STIX]{x1D70F}$ via

with

Defining

with $\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ and $(v_{\Vert }\boldsymbol{n}+\boldsymbol{v}_{d})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ slow functions of space, gives upon integrating along the trajectory

where the weak $\unicode[STIX]{x1D6E9}(\unicode[STIX]{x1D70F})$ and/or $\unicode[STIX]{x1D70F}$ dependence of the coefficient $W_{p,m}\exp (\text{i}p\unicode[STIX]{x1D6FC})$ will turn out to be unimportant because of our approximation in the next paragraph. Integrating from $\unicode[STIX]{x1D70F}\rightarrow -\infty$ to $\unicode[STIX]{x1D70F}=0$ , and assuming any memory of the initial state is lost so we can take $g_{p}[\unicode[STIX]{x1D6E9}(\unicode[STIX]{x1D70F}\rightarrow -\infty )]=0$ gives the general form

The usual derivations of the QL operator always neglect any memory of the initial condition of the charges. The use of unperturbed trajectories might seem to imply a long memory of the past, in which case even collisions might matter. However, in practice the applied rf fields actually perturb these trajectories and thereby introduce stochasticity (Becoulet et al. Reference Becoulet, Gambier and Samain1991). As a result, for a QL treatment the unperturbed trajectories are only approximately valid between successive passes through resonance.

To simplify this expression we integrate by parts by taking advantage of the rapid variation of the exponential compared to the slow variation of $W_{p,m}$ . To do so we use

The

term includes the following contributions:

and

with $P$ the total plasma pressure.

We could also integrate (4.63) over the periodic poloidal motion of the trapped and passing particles by assuming that the unperturbed trajectories remained valid for many complete poloidal periods. This procedure could imply weaker stochasticity since successive passes through the resonance are correlated, leading to bounce and poloidal transit resonances appearing in the resonant denominator. Such behaviour is inherent in transit averaged descriptions which make an implicit assumption of weaker stochasticity.

The coefficient $W_{p,m}$ is expected to be slowly varying compared to the variation in the vicinity of the denominator since $W_{p,m}$ is constant in a constant magnetic field. Assuming the denominator changes slowly as the charged particle passes through resonance,

then an additional integration by parts yields the form convenient for full wave treatments:

where we neglect the small corrections. For unmagnetized plasmas or low frequency waves ( $\unicode[STIX]{x1D714}\ll \unicode[STIX]{x1D6FA}$ ) only the $p=0$ term matters since all other terms will be small by $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FA}$ . For higher frequencies such that $\unicode[STIX]{x1D714}\sim p\unicode[STIX]{x1D6FA}+k_{\Vert }v_{\Vert }+\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ only the $p$ term in the sum is significant since all other terms will be small by $|\unicode[STIX]{x1D714}-p\unicode[STIX]{x1D6FA}-k_{\Vert }v_{\Vert }-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}|/p\unicode[STIX]{x1D6FA}\ll 1$ . For higher frequencies, the $k_{\Vert }v_{\Vert }$ and $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ terms are important because they determine the location of the wave damping – and this can be quite different than where $\unicode[STIX]{x1D714}=p\unicode[STIX]{x1D6FA}$ . This form is used to evaluate the density and current in full wave treatment since we can write $g$ in drift kinetic variables as

This form is not a strict Fourier representation in poloidal angle since the coefficients depend on poloidal angle. To see that this is a compact and efficient means of treating the poloidal variation of the resonant part of the distribution function, assume that there is only a single poloidal mode for the applied rf electric field of (4.21). Then there is only a single $m$ term in our non-standard representation for the distribution function $g$ of (4.71). Due to the resonant denominator its coefficient varies strongly with poloidal angle $\unicode[STIX]{x1D717}$ . As a result, a standard Fourier representation would require an infinite sum of poloidal modes.

The preceding form (4.71) is convenient to use when forming the density and species current,

and

with

and

where Re denotes the real part. As we did not Fourier decompose $g_{p}$ in poloidal angle, our full wave results are quite general and compact, with all poloidal angle coupling due to magnetic field variation retained. We have allowed

so that drift effects can become significant for the ions (denoted by a subscript $i$ ) when $k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$ for $B_{p}\ll B$ . Here $\unicode[STIX]{x1D70C}_{pi}=\unicode[STIX]{x1D70C}_{i}B/B_{p}$ is the ion poloidal gyroradius and $B_{p}$ the poloidal magnetic field. These moment expressions ignore any drift corrections from $\unicode[STIX]{x1D6FF}_{p,m}$ implying that $\unicode[STIX]{x1D714}\gg \unicode[STIX]{x1D6E5}_{m}$ for $p=0$ and $\unicode[STIX]{x1D714}-k_{\Vert }v_{\Vert }-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}\sim p\unicode[STIX]{x1D6FA}\gg \unicode[STIX]{x1D6E5}_{m}$ for $p\neq 0$ .

To form the quasilinear (QL) operator, the principle value part of the integral

is neglected. The preceding expression uses the $\unicode[STIX]{x1D6FF}$ function form found from a causal Laplace transform with a positive imaginary part to $\unicode[STIX]{x1D714}$ , $\text{Im}\unicode[STIX]{x1D714}>0$ ,

Using this result gives the approximation needed for the QL operator to be

For Kennel & Engelmann (Reference Kennel and Engelmann1966) to be valid we need to assume $k_{\Vert }R\gg k_{\bot }\unicode[STIX]{x1D70C}$ .

To examine the validity of the preceding expressions we use

Then the width of the resonance is determined by the broadening so we can estimate

As a result, we expect our results to be valid when

since poloidal angle variation of $W_{p,m}$ varies slowly through the narrow resonance region that is much less than the minor radius $r$ appearing in the coefficient of the poloidal varying portion of the magnetic field $B\propto R^{-1}\propto [1+(r/R)\cos \unicode[STIX]{x1D717}]^{-1}$ . However, trapped particles reflecting ( $v_{\Vert }=0$ ) in the resonance layer $\unicode[STIX]{x1D714}\simeq p\unicode[STIX]{x1D6FA}$ are known to violate (4.82) when $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}\rightarrow 0$ (Catto & Myra Reference Catto and Myra1992; Belikov & Kolesnichenko Reference Belikov and Kolesnichenko1994; Catto, Myra & Russell Reference Catto, Myra and Russell1994). Whether this continues to be the case when $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ is retained is less clear since $\boldsymbol{v}_{d}\neq 0$ when $v_{\Vert }=0$ .

The new term we obtain in the resonant denominator or argument of the delta function is the expected Doppler broadening of the resonance due to the shift $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ caused by magnetic drift effects. As a result of this drift motion the actual resonance is shifted to a different location on the flux surface with the resonance occurring at a slightly different poloidal location in velocity space experiencing different wave fields. The distinction between the gyrofrequency of the guiding centre and charge does not matter since only $k_{\bot }\unicode[STIX]{x1D70C}_{i}\gg 1$ are of interest. The damping associated with the $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ magnetic drift is illustrated in Bajaj & Krall (Reference Bajaj and Krall1972) and Lee et al. (Reference Lee, Myra and Catto1983) for ion drift cyclotron waves. Although unimportant for our purposes, there is also Landau damping associated with the distinction between the guiding centre and charge location as illustrated in Antonsen & Manheimer (Reference Antonsen and Manheimer1978) and Lee et al. (Reference Lee, Myra and Catto1983) for ion Bernstein modes. All these new broadening terms are expected to matter more for ions than electrons.

5 Quasilinear operator

To form the QL operator we need $g$ in drift kinetic variables. Expanding the exponential we obtain

Then, expanding the exponential in (4.79)

gives

where we have carefully kept the distinction between drift kinetic and gyrokinetic variables in the exponential to get the $L$ factor. Consistent with our orderings we ignore the distinction between gyrokinetic and drift kinetic variables in the gyrofrequency.

In the QL operator we can use the delta function to simplify $W_{p,m}$ by letting

Moreover, we need to account for the difference between $f_{1}$ and $g$ ,

The difference $f_{1}-g$ is gyrophase independent, while $\langle \boldsymbol{v}_{\bot }\rangle _{\unicode[STIX]{x1D711}}=0$ , leaving

where we use $\unicode[STIX]{x1D735}_{v}\unicode[STIX]{x1D711}=v_{\bot }^{-2}\boldsymbol{n}\times \boldsymbol{v}$ to obtain the lowest-order result $\unicode[STIX]{x1D735}_{v}v\times \unicode[STIX]{x1D735}_{v}\unicode[STIX]{x1D707}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}\unicode[STIX]{x1D711}=1/vB$ .

We must again consider $\unicode[STIX]{x1D707}$ to higher order to retain drift effects

and

giving

To simplify further for harmonic fields $\unicode[STIX]{x1D73A}$ , $\unicode[STIX]{x1D74C}$ and $\unicode[STIX]{x1D747}$ we use

and

so if $\unicode[STIX]{x1D747}_{m}=\boldsymbol{e}_{m}$ and $\unicode[STIX]{x1D74C}_{m}=\text{i}\boldsymbol{e}_{m}$ we get zero. The $\langle \ldots \rangle _{\unicode[STIX]{x1D717}}$ in $\langle \ldots \rangle _{cg}$ replaces the double poloidal angle sum over $m$ and $m^{\prime }$ by a single sum over $m$ by introducing the Kronecker delta function $\unicode[STIX]{x1D6FF}_{m^{\prime }m}$ . By using (5.10) and (5.11) we are not allowing poloidal mode numbers that differ by order unity to contribute to the quasilinear operator. This property of $\langle \ldots \rangle _{\unicode[STIX]{x1D717}}$ is necessary to maintain the correct entropy production for our QL operator and recover the Kennel–Engelmann (Reference Kennel and Engelmann1966) operator. The preceding averages are performed in real physical space, unlike the elegant Hamiltonian action-angle treatment of Kaufman (Reference Kaufman1972) and Eriksson & Helander (Reference Eriksson and Helander1994) that employ similar averages in phase space. Their Hamiltonian treatments lead to transit average quasilinear operators in action-angle variables, while we will find the non-transit averaged form in configuration space variables.

Recalling that for $\text{Im}\,\unicode[STIX]{x1D714}>0$

we find

Therefore, we need only consider

Using the preceding and

and

gives

In addition, we use

with

where in the last form we have neglected unimportant drift corrections. Using this result along with earlier results, and $\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{e}_{\bot p}v_{\bot }=p\unicode[STIX]{x1D6FA}\text{J}_{p}(w)$ and $\unicode[STIX]{x1D714}-p\unicode[STIX]{x1D6FA}-k_{\Vert }v_{\Vert }-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}=0$ gives

Then we need only note that all the terms in

are small compared to the terms just evaluated.

Combining the preceding results, we find our final form for the QL operator to be

with

manifestly positive. Consequently, we have managed to carefully derive the desired tokamak form for the QL operator that is the same as the Kennel–Englemann (Reference Kennel and Engelmann1966) form except for the simple $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{d}$ modification to the argument of the delta function. Our QL ordering $Q\sim \unicode[STIX]{x1D708}f_{0}$ with $D\sim Z^{2}e^{2}|\boldsymbol{e}_{m}|^{2}/M^{2}\unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D714}v_{\text{quiv}}^{2}$ and $v_{\text{quiv}}\sim Ze|\boldsymbol{e}_{m}|/M\unicode[STIX]{x1D714}$ , then requires small quiver speeds $v_{\text{quiv}}$ compared to the thermal speed $v_{\text{thermal}}$ such that $v_{\text{quiv}}^{2}/v_{\text{thermal}}^{2}\sim \unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}\ll 1$ .

The preceding form satisfies a QL entropy principle. Multiplying Q by $\ell nf_{0}$ and integrating over all velocity space using $\text{d}^{3}v\rightarrow 2\unicode[STIX]{x03C0}vB\,\text{d}v\,\text{d}\unicode[STIX]{x1D707}/|v_{\Vert }|$ gives the entropy production to be negative definite as required,

By changing to the cylindrical velocity space variables $v_{\bot }$ , $v_{\Vert }$ , and $\unicode[STIX]{x1D711}$ we can rewrite the quasilinear operator in vector form as

with

where the various diffusivities are defined as

and

For simulations the delta function can be approximated as the limiting form of a Gaussian. Also, using form (5.25), we can easily evaluate the energy moment of $Q\{\,f_{0}\}$ as

For completeness we also note that the transit averaged QL operator is

with $\overline{D}$ the transit average of the full diffusivity

Sometimes a pitch angle variable is used in place of $\unicode[STIX]{x1D707}$ . The transit average of $D$ is most easily performed explicitly by using the argument of the delta function to select the poloidal angle at resonance.

Our expression for the transit averaged quasilinear operator $\overline{Q\{\,f_{0}\}}$ is equivalent to the usual Kennel & Engelmann (Reference Kennel and Engelmann1966) and Eriksson & Helander (Reference Eriksson and Helander1994) result written in an explicit and compact form in drift kinetic variables, but without the radial diffusion terms. It is consistent with the transit averaged form of Belikov & Kolesnichenko (Reference Belikov and Kolesnichenko1994) that extended the result of Catto & Myra (Reference Catto and Myra1992) to allow finite $k_{\Vert }$ and retain radial diffusion. Both references explicitly treat the distinction between trapped and passing particles, but without the drift modifications considered here.

In heating and current drive applications the wave fields for the QL diffusivities are assumed to be known from full wave codes such as TORIC (Brambilla Reference Brambilla1999; Wright et al. Reference Wright, Bonoli, Brambilla, Meo, Dazevedo, Batchelor, Jaeger, Berry, Phillips and Pletzer2004) and AORSA (Jaeger et al. Reference Jaeger, Berry, Dazevedo, Batchelor, Carter and White2002). The density and currents needed for a full wave code have already been partially evaluated by performing the required gyroaverages in (4.72) to (4.75). The remaining integrals can be performed once $f_{0}$ is found numerically by solving the QL equation. Coarse grain averaging is unnecessary for the full wave representation of the linear wave fields. Otherwise they are derived with the same approximations as used to obtain the QL operator and so provide a consistent evaluation for a full wave code in a representation that is expected to require keeping fewer poloidal modes. They can be converted to $v_{\Vert }$ and $\unicode[STIX]{x1D707}$ variables so the non-transit averaged form of the kinetic equation for $f_{0}$ can be employed to keep streaming effects that can sometimes play a role (Helander & Catto Reference Helander and Catto2001).

6 Summary

We have derived the modification of the standard Kennel–Engelmann (Reference Kennel and Engelmann1966) QL operator that incorporates magnetic drift effects in tokamak geometry in a non-transit averaged form and continues to have non-negative entropy production. The only modification is in the argument of delta function that now includes a broadening term due to the total magnetic drift as well as usual Landau and cyclotron broadening of the resonance. To carefully derive this drift modified form for the quasilinear operator is not as straightforward as might be expected because the higher-order magnetic moment must be retained. The need to go to higher order requires the unconventional steps of solving explicitly for two non-resonant terms in the perturbed distribution function that do not alter the QL operator and directly lead to its most compact form. Our QL form also has the advantages that it is explicitly in drift kinetic variables and not transit averaged.

The perturbed distribution function need not be periodic in poloidal angle because of the presence of trapped particles. Therefore, it is not Fourier decomposed in poloidal angle like the periodic wave fields. Consequently, a compact modified Fourier representation for the perturbed distribution function is found having poloidally angle dependent coefficients. This form is expected to be useful for full wave treatments.