Electromagnetic zonal flow residual responses

TL;DR

This paper extends the calculation of electromagnetic zonal flow residuals in tokamak plasmas, providing new formulas and insights into the behavior of electromagnetic perturbations and their relaxation properties.

Contribution

It generalizes the collisionless zonal flow residual calculation to include electromagnetic effects with fully self-consistent fields, offering new test expressions for gyrokinetic codes.

Findings

Derived simple expressions for electromagnetic residual responses.

Showed that certain magnetic perturbations do not relax to flux functions.

Provided a comprehensive electromagnetic test for gyrokinetic simulations.

Abstract

The collisionless axisymmetric zonal flow residual calculation for a tokamak plasma is generalized to include electromagnetic perturbations. We formulate and solve the complete initial value zonal flow problem by retaining the fully self-consistent axisymmetric spatial perturbations in the electric and magnetic fields. Simple expressions for the electrostatic, shear and compressional magnetic residual responses are derived that provide a fully electromagnetic test of the zonal flow residual in gyrokinetic codes. Unlike the electrostatic potential, the parallel vector potential and the parallel magnetic field perturbations need not relax to flux functions for all possible initial conditions.

Equations

$$\mathbf{B}={\mathbf{B}}_0+{\mathbf{B}}_1,$$

$${\mathbf{B}}_0=I(\psi )\mathrm{NewA}\zeta +\mathrm{NewA}\zeta \times \mathrm{NewA}\psi =B_0\mathbf{b}$$

$${\mathbf{B}}_1=\mathrm{NewA}\times \mathbf{A},$$

$$\mathbf{A}=A_{\Vert }\mathbf{b}+A_{\perp }{\mathbf{e}}_{\perp }+A_{\times }{\mathbf{e}}_{\times }.$$

$$\left\{A_{\Vert },A_{\perp },A_{\times }\right\}=\left\{{\tilde{A}}_{\Vert },{\tilde{A}}_{\perp },{\tilde{A}}_{\times }\right\}\exp [iS(\psi )],$$

$${\mathbf{B}}_1\simeq B_{\Vert }\mathbf{b}-i{k}_{\perp }{A}_{\Vert }{\mathbf{e}}_{\times },$$

$${\mathbf{B}}_1=\delta \mathrm{NewA}\zeta -\mathrm{NewA}\zeta \times \mathrm{NewA}A,$$

$$\mathrm{NewA}A\simeq \frac{i{k}_{\perp }A}{R{B}_p}\nabla \psi ,$$

$${\mathbf{B}}_1\simeq \delta \mathrm{NewA}\zeta -\frac{i{k}_{\perp }A}{R{B}_p}\mathrm{NewA}\zeta \times \mathrm{NewA}\psi .$$

$${\mathbf{e}}_{\times }\cdot \mathrm{NewA}\zeta =-\frac{B_p}{R{B}_0}$$

$${\mathbf{e}}_{\times }\cdot (\mathrm{NewA}\zeta \times \mathrm{NewA}\psi )=\frac{I{B}_p}{R{B}_0},$$

$$B_{\Vert }={\mathbf{B}}_1\cdot \mathbf{b}=\frac{I}{R^2{B}_0}\delta -\frac{i{k}_{\perp }{B}_p}{R{B}_0}A$$

$$A_{\Vert }=-\frac{1}{i{k}_{\perp }}{\mathbf{B}}_1\cdot {\mathbf{e}}_{\times }=\frac{B_p}{i{k}_{\perp }R{B}_0}\delta +\frac{I}{R^2{B}_0}A.$$

$$\delta =R^2{\mathbf{B}}_1\cdot \mathrm{NewA}\zeta =\frac{I}{B_0}{B}_{\Vert }+\frac{i{k}_{\perp }R{B}_p}{B_0}{A}_{\Vert }$$

$$A=-\frac{R}{i{k}_{\perp }{B}_p}{\mathbf{B}}_1\cdot (\mathrm{NewA}\zeta \times \mathrm{NewA}\psi )=-\frac{R{B}_p}{i{k}_{\perp }{B}_0}{B}_{\Vert }+\frac{I}{B_0}{A}_{\Vert }.$$

$$q=q_0+q_1,$$

$$q_0=\frac{1}{2\pi }{\int }_0^{2\pi }\frac{{\mathbf{B}}_0\cdot \mathrm{NewA}\zeta }{{\mathbf{B}}_0\cdot \mathrm{NewA}\theta }d\theta$$

$$q_1=\frac{1}{2\pi }{\int }_0^{2\pi }\left(\frac{{\mathbf{B}}_1\cdot \mathrm{NewA}\zeta }{{\mathbf{B}}_0\cdot \mathrm{NewA}\zeta }-\frac{{\mathbf{B}}_1\cdot \mathrm{NewA}\theta }{{\mathbf{B}}_0\cdot \mathrm{NewA}\theta }\right)\frac{{\mathbf{B}}_0\cdot \mathrm{NewA}\zeta }{{\mathbf{B}}_0\cdot \mathrm{NewA}\theta }d\theta .$$

$$q_1=\frac{1}{2\pi }\frac{i{k}_{\perp }}{R{B}_p}{\int }_0^{2\pi }\frac{A_{\Vert }}{\mathbf{b}\cdot \mathrm{NewA}\theta }d\theta =\frac{1}{2\pi }\frac{i{k}_{\perp }}{R{B}_p}{\int }_0^{2\pi }\left(\frac{B_p\delta }{i{k}_{\perp }R}+\frac{IA}{R^2}\right)\frac{d\theta }{{\mathbf{B}}_0\cdot \mathrm{NewA}\theta }.$$

$$f_0=n{\left(\frac{M}{2\pi T}\right)}^{3/2}\exp \left(-\frac{M{v}^2}{2T}\right),$$

$${\dot{f}}_1=\frac{\partial f_1}{\partial t}+\mathbf{v}\cdot \mathrm{NewA}{f}_1+\Omega \mathbf{v}\times \mathbf{b}\cdot {\nabla }_v{f}_1=\frac{Ze}{T}{\mathbf{E}}_1\cdot \mathbf{v}{f}_0,$$

$$\dot{h}=\frac{\partial h}{\partial t}+\mathbf{v}\cdot \mathrm{NewA}h+\Omega \mathbf{v}\times \mathbf{b}\cdot {\nabla }_{v}h=\frac{Ze}{T}\left[\frac{\partial \Phi }{\partial t}-\frac{\partial }{\partial t}\left(\frac{\mathbf{v}\cdot \mathbf{A}}{c}\right)\right]f_0.$$

$${\psi }_{\ast }=\psi -\frac{Mc}{Ze}R\hat{\zeta }\cdot \mathbf{v}=\psi +\frac{1}{\Omega }{\mathbf{v}}_{\perp }\times \mathbf{b}\cdot \mathrm{NewA}\psi -\frac{I{v}_{\Vert }}{\Omega }.$$

$$\Phi =\tilde{\Phi }\exp [iS(\psi )],$$

$$\frac{\partial h}{\partial t}+v_{\Vert }\mathbf{b}\cdot \mathrm{NewA}h=\frac{Ze{f}_0}{cT}{⟨\exp [iS(\psi )]\frac{\partial }{\partial t}\left[c\tilde{\Phi }-v_{\Vert }{\tilde{A}}_{\Vert }-{\mathbf{v}}_{\perp }\cdot \left({\mathbf{e}}_{\perp }{\tilde{A}}_{\perp }+{\mathbf{e}}_{\times }{\tilde{A}}_{\times }\right)\right]⟩}_{\phi },$$

$$h=h({\psi }_{\ast },\theta ,v,\mu ,t)={\tilde{h}}_{\ast }({\psi }_{\ast },\theta ,v,\mu ,t)\exp \left[iS({\psi }_{\ast })\right],$$

$$\frac{\partial {\tilde{h}}_{\ast }}{\partial t}+v_{\Vert }\mathbf{b}\cdot \mathrm{NewA}{\tilde{h}}_{\ast }=\frac{Ze{f}_0}{cT}{⟨\exp [iQ-iL]\frac{\partial }{\partial t}\left[c\tilde{\Phi }-v_{\Vert }{\tilde{A}}_{\Vert }-{\mathbf{v}}_{\perp }\cdot \left({\mathbf{e}}_{\perp }{\tilde{A}}_{\perp }+{\mathbf{e}}_{\times }{\tilde{A}}_{\times }\right)\right]⟩}_{\phi }.$$

$$\frac{\partial {\tilde{h}}_{\ast }}{\partial t}+v_{\Vert }\mathbf{b}\cdot \mathrm{NewA}{\tilde{h}}_{\ast }=\frac{Ze{f}_0}{cT}\exp (iQ)\frac{\partial }{\partial t}\left[J_0\left(\frac{k_{\perp }{v}_{\perp }}{\Omega }\right)(c\tilde{\Phi }-v_{\Vert }{\tilde{A}}_{\Vert })+\frac{2J_1(k_{\perp }{v}_{\perp }/\Omega )}{k_{\perp }{v}_{\perp }/\Omega }\frac{v_{\perp }^2}{2\Omega }{\tilde{B}}_{\Vert }\right],$$

$$\frac{\partial \overline{{\tilde{h}}_{\ast }}}{\partial t}=\frac{Ze{f}_0}{cT}\frac{\partial }{\partial t}\overline{\left[J_0\left(\frac{k_{\perp }{v}_{\perp }}{\Omega }\right)(c\tilde{\Phi }-v_{\Vert }{\tilde{A}}_{\Vert })+\frac{2J_1(k_{\perp }{v}_{\perp }/\Omega )}{k_{\perp }{v}_{\perp }/\Omega }\frac{v_{\perp }^2}{2\Omega }{\tilde{B}}_{\Vert }\right]\exp (iQ)},$$

$$\overline{{\tilde{h}}_{\ast }}=\overline{{\tilde{h}}_{\ast }(t=0)}+$$

Full text

††volume: ?

{PACS}

?

1 Introduction

A zonal flow is a sheared flow generated by turbulence that has small scale structure compared to the system size in the radial direction and is global in extent in the other directions. In a tokamak an electrostatic zonal flow appears as a radially varying electric field drift due to a radial electric field with rapid radial variation, but with no toroidal variation. It helps reduce and regulate the turbulent transport level in tokamaks through shear-enhanced decorrelation of turbulent structures (Biglari et al., 1990; Terry, 2000). Its importance was discovered when a discrepancy between gyrokinetic and gyrofluid descriptions of ion temperature gradient (ITG) turbulence was observed in the earliest nonlinear, electrostatic, $\delta f$, flux-tube, particle-in-cell code (now called PG3EQ) (Dimits et al., 1996), where $\delta f$ is the perturbation way from the Maxwellian. The key role of zonal flow in controlling and reducing ITG turbulent transport, especially near marginal stability, soon became apparent. Insights into zonal flow behavior (missed in early gyrofluid codes) came from code simulations and comparisons (Dimits et al., 2000), leading to an understanding that there was a nonlinear Dimits shift away from the ITG linear stability threshold (Dimits et al., 1996). Rosenbluth & Hinton (1998) developed an electrostatic analytic check to show that the zonal flow damps to a non-zero residual level in a collisionless, axisymmetric plasma due to polarization effects associated with the magnetic drifts.

The standard zonal flow residual calculations are electrostatic and assume axisymmetry is maintained. An initial value problem is solved to find the residual zonal flow level once any initial poloidal angle dependence in the electrostatic potential is temporally damped away via the geodesic acoustic mode (GAM) (Winsor et al., 1968). The initial condition must normally be chosen to depend on poloidal angle to generate a transiently evolving zonal flow or else it will be a homogeneous solution to the non-transit averaged drift kinetic or gyrokinetic equation. The initial distribution function is not allowed to depend on gyro-phase since it is assumed that any initial transient response associated with the fast gyro-motion has already damped to its classical polarization level. The residual zonal flow level has proven to be an important electrostatic test of gyrokinetic codes in general (Dimits et al., 2000) and the GS2 code in particular (Xiao et al., 2007a), including collisional damping (Hinton & Rosenbluth, 1999; Xiao et al., 2007b), short wavelength effects (Jenko et al., 2000; Xiao & Catto, 2006b) and the effect of shaping (Belli, 2006; Xiao & Catto, 2006a).

In subsequent sections we generalize this axisymmetric electrostatic model to its fully electromagnetic counterpart for a tokamak. Unlike the procedure of Terry et al. (2013), which treats the effect of an externally imposed, stationary (non-evolving) and non-axisymmetric radial magnetic field perturbation on an equilibrium, we formulate and solve a description retaining the fully self-consistent axisymmetric spatial perturbations in the magnetic field. The poloidal dependence of the perturbed shear and compressional magnetic field perturbations are retained as drives in the kinetic equation, quasineutrality, and the parallel and perpendicular components of Ampère’s law. The system of equations are then solved to obtain the complete self-consistent response within a Vlasov-Maxwell description. The expressions obtained by solving this initial value problem provide 15 fully electromagnetic tests of the zonal flow residual in gyrokinetic codes. The poloidally dependent initial conditions for the fields and the distribution function are chosen to satisfy quasineutrality and Ampère’s law at $t=0$. We assume any GAM behaviour due to poloidal variation has damped away so that only the residual zonal flow levels are obtained. Importantly, the residual zonal flow levels must allow for poloidal variation of the parallel vector potential and parallel magnetic field perturbations for all initial conditions. The description is general enough that even in the absence of any initial electrostatic perturbation, a magnetic perturbation is able to generate a zonal flow response.

The subsequent sections are organized as follows. First, in section 2 we specify the representations of the perturbed and unperturbed fields. The kinetic equation is given in section 3, then a suitable initial condition in terms of perturbed fields and distribution function is chosen in section 4. The system is closed with Maxwell’s equations and solved in section 5. Approximate expressions for the zonal flow responses in the various fields are given in section 6, before we briefly summarize our results in section 7.

2 Potentials, fields, and currents

The standard electrostatic zonal flow residual calculation (Rosenbluth & Hinton, 1998) assumes axisymmetry is preserved during the time evolution of the zonal flow. We seek to generalize this model to its fully electromagnetic counterpart in a tokamak, assuming that the magnetic field remains axisymmetric at all times. Unperturbed quantities are assumed to evolve slowly compared to the zonal flow relaxation.

The total magnetic field is

[TABLE]

where

[TABLE]

is the background axisymmetric magnetic field, and $\mathbf{B}_{1}$ is the perturbed magnetic field. Here $I=I(\psi)$ must be a flux function to make the unperturbed radial current density vanish, $2\pi\psi$ is the unperturbed poloidal flux, $R{\rm NewA}\zeta=\hat{\zeta}$ is the toroidal unit vector, with $R$ the major radius, $B_{0}=|\mathbf{B}_{0}|$ is the magnitude of the unperturbed magnetic field, and $\mathbf{b}=\mathbf{B}_{0}/B_{0}$ is the unit vector in the direction of the unperturbed magnetic field.

We start by representing $\mathbf{B}_{1}$ in a form convenient to derive the gyroaveraged kinetic equation,

[TABLE]

where

[TABLE]

Here $\mathbf{e}_{\perp}={\rm NewA}\psi/RB_{p}$, $\mathbf{e}_{\times}=\mathbf{b}\times\mathbf{e}_{\perp}$, $|{\rm NewA}\psi|=RB_{p}$ and $B_{p}$ is the unperturbed poloidal magnetic field. We assume that the characteristic length scale of $\mathbf{B}_{1}$ perpendicular to the background magnetic field $\mathbf{B}_{0}$ is small compared to the characteristic size of the device. Since we are considering zonal components, the perpendicular gradient is mainly in the radial direction. To describe this rapid radial variation of the perturbed fields, we use the eikonal form

[TABLE]

with ${\rm NewA}S=S^{\prime}{\rm NewA}\psi=k_{\perp}\mathbf{e}_{\perp}$. The coefficients with tilde are functions of time and are only allowed to be slow functions of $\psi$ and $\theta$, for example, varying as $\cos\theta$. Using equations (3), (4) and (5), we obtain

[TABLE]

where $B_{\|}=ik_{\perp}A_{\times}$ is the parallel component of $\mathbf{B}_{1}$. The form for $\mathbf{B}_{1}$ in (6) is the most common way to express the perturbed magnetic field in gyrokinetic simulations. Note that the component $A_{\perp}$ of the vector potential never appears in the final expression for $\mathbf{B}_{1}$, and can be safely ignored.

The form for $\mathbf{B}_{1}$ in (6) ensures that the magnetic field $\mathbf{B}_{1}$ is axisymmetric. We can make this more explicit by showing that (6) is equivalent to the axisymmetric form

[TABLE]

where $\delta(\psi,\theta)$, which need not be a flux function, is the perturbation to $I(\psi)$, and $-A(\psi,\theta)$ is the perturbation to $\psi$. To match equations (6) and (7), we must realize that $A$ has an eikonal form similar to those in (5). Thus,

[TABLE]

and as a result, equation (7) gives

[TABLE]

Equation (9) proves that $\mathbf{B}_{1}$ is, to lowest order in $(k_{\perp}a)^{-1}\ll 1$, parallel to the flux surface. Then, it is easy to obtain the relations between $\delta$, $A$, $B_{\|}$ and $A_{\|}$. From (6), we deduce that $B_{\|}=\mathbf{B}_{1}\cdot\mathbf{b}$ and $A_{\|}=-\mathbf{B}_{1}\cdot\mathbf{e}_{\times}/ik_{\perp}$. Substituting into these two equations the form for $\mathbf{B}_{1}$ given in (9), and using

[TABLE]

and

[TABLE]

we obtain

[TABLE]

and

[TABLE]

Similarly, from (9), we find $\delta=R^{2}\mathbf{B}_{1}\cdot{\rm NewA}\zeta$ and $A=-(R/ik_{\perp}B_{p})\mathbf{B}_{1}\cdot({\rm NewA}\zeta\times{\rm NewA}\psi)$. Substituting into these equations the form for $\mathbf{B}_{1}$ given in (6), and using (10) and (11), we obtain

[TABLE]

and

[TABLE]

Expressions (12), (13), (14) and (15) allow us to change from the form most convenient for gyrokinetics, in (6), to the axisymmetric form in (9).

It is of interest to discuss the possible changes that the magnetic field can undergo. It is possible to change the direction of the magnetic field lines without changing the magnitude of $\mathbf{B}$ if $B_{\|}=0=I\delta-ik_{\perp}RB_{p}A$. To avoid changing the direction of the field lines, the perturbation must satisfy $\mathbf{B}_{1}\times\mathbf{B}_{0}=0$, or $A_{\|}=0=RB_{p}\delta+ik_{\perp}IA$. Changing the local theta dependent direction of the magnetic field line does not necessarily imply a change in the safety factor $q(\psi)=(2\pi)^{-1}\int_{0}^{2\pi}(\mathbf{B}\cdot{\rm NewA}\zeta/\mathbf{B}\cdot{\rm NewA}\theta)d\theta$. To verify this we write the safety factor as

[TABLE]

where

[TABLE]

is the unperturbed safety factor and

[TABLE]

is the result of the perturbation. Note that our poloidal angle-like variable $\theta$ is not changed by the perturbations. Using (6) and (9), and recalling that $k_{\perp}/RB_{p}$ is independent of $\theta$, the perturbation $q_{1}$ becomes

[TABLE]

When $q_{1}$ changes, the lines in the flux surface change topology by switching between rational and irrational – the only form of reconnection allowed for axisymmetric perturbations.

3 Kinetic equation and solution

We need to solve the linearized gyrokinetic equation in which the unperturbed quantities are time independent or evolve slowly compared to the zonal flow relaxation. The unperturbed ion distribution function $f_{0}$ is assumed to be Maxwellian:

[TABLE]

with $n$, $M$ and $T$ the ion density, mass and temperature, respectively. Then, as we shall consider a collisionless plasma, the linearized distribution function $f_{1}$ satisfies the Vlasov equation

[TABLE]

where $\Omega=ZeB_{0}/Mc$ with $Ze$ the ion charge and $c$ the speed of light. We have neglected the unperturbed electric field. For the perturbed electric field we use $\mathbf{E}_{1}=-{\rm NewA}\Phi-c^{-1}\partial\mathbf{A}/\partial t$, with $\Phi$ the perturbed electrostatic potential. To remove the adiabatic piece we let $f_{1}=h-(Ze\Phi/T)f_{0}$ to obtain

[TABLE]

This is the form that we will use to obtain the desired gyrokinetic equation.

Rather than perform a conventional gyrokinetic treatment (Catto, 1978) of (22), we use the canonical angular momentum $\psi_{\ast}$ for our radial variable when we change variables (Kagan & Catto, 2008). Then using $\mathbf{v}=\mathbf{v}_{\perp}+v_{\|}\mathbf{b}$ gives

[TABLE]

As for the vector potential in (5), we describe the rapid radial variation of the perturbed electrostatic potential $\Phi$ by the eikonal expression

[TABLE]

where the coefficient $\tilde{\Phi}$ is a function of time and is only allowed to be a slow functions of $\theta$ and $\psi$. Changing from $\psi$, $\theta$, $\zeta$ and $\mathbf{v}$ variables to $\psi_{\ast}$, $\theta$, $\zeta$, $v$, $\mu=v_{\perp}^{2}/2B_{0}$, and gyro-phase $\varphi$, with $\mathbf{v}_{\perp}=v_{\perp}(\mathbf{e}_{\perp}\cos\varphi+\mathbf{e}_{\times}\sin\varphi)$, and using $\dot{\psi}_{\ast}=0$ to remove the $\psi_{\ast}$ derivative (Kagan & Catto, 2009) yields the lowest order gyrokinetic equation

[TABLE]

where the $\theta$ dependence of $h$ is assumed slow, drift corrections to parallel streaming are neglected as small, and the gyroaverage $\langle\dots\rangle_{\varphi}=(2\pi)^{-1}\oint d\varphi(\dots)$ is performed at fixed $\psi_{\ast}$.

Next we write $h$ in the eikonal form

[TABLE]

where only a weak $\psi_{\ast}$ dependence of $\tilde{h}_{\ast}$ not captured by $S(\psi_{\ast})$ is allowed, $v=|\mathbf{v}|$, and $\mu=v_{\perp}^{2}/(2B_{0})$. Then, we Taylor expand to obtain $S(\psi)-S(\psi_{\ast})=Q-L+\dots$, where $L=(k_{\perp}v_{\perp}/\Omega)\sin\varphi$, $Q=S^{\prime}Iv_{\|}/\Omega=k_{\perp}v_{\|}/\Omega_{p}$, and $\Omega_{p}=ZeB_{p}/Mc=\Omega B_{p}/B_{0}$. We retain the order $\epsilon$ Shafranov shift $\Delta_{S}$ of the flux surfaces by writing $R=R_{0}(\psi)+r(\psi)\cos\theta$ with $R_{0}(\psi=0)=R_{0}(0)$ the location of the magnetic axis, $R_{0}(\psi)=R_{0}(0)-\Delta_{S}$, and $r$ the minor radius for circular flux surfaces, then $\epsilon=r/R_{0}(\psi)$ is a flux function and the ratio of the poloidal over the toroidal magnetic field is $B_{p}/B_{t}=(\epsilon/q)[1-\Delta_{S}\cos\theta+\mathcal{O}(\epsilon^{2})]$ with $I=RB_{t}$. As a result, the lowest order gyrokinetic equation becomes

[TABLE]

Performing the gyroaverages we obtain the desired form of the ion gyrokinetic equation

[TABLE]

where $\tilde{B}_{\|}=ik_{\perp}\tilde{A}_{\times}$, $\mathbf{b}\cdot{\rm NewA}=(qR)^{-1}\partial/\partial\theta$, the coefficient of the $\tilde{A}_{\perp}$ term has gyroaveraged to zero, and $J_{0}$ and $J_{1}$ denote Bessel functions of the first kind.

To lowest order the streaming term dominates for a weakly collisional plasma. The damping away of any initial $\theta$ dependence leads to the GAM behaviour observed during the early evolution to the final residual zonal flow steady state. We are not interested in the GAM behaviour so we annihilate the $\theta$ derivative in (28) to obtain the transit averaged gyrokinetic equation

[TABLE]

where $\overline{\tilde{h}_{\ast}}$ is the $\theta$ independent long time solution, while $\tilde{h}_{\ast}(t=0)$ is allowed to depend on $\theta$ so we can relate an initial $\theta$ dependent perturbation to the final steady state $\theta$ independent solution. The transit average of any quantity $X$ is defined as $\overline{X}=\oint d\tau X/\oint d\tau$, with $d\tau=d\theta/(v_{\|}\mathbf{b}\cdot{\rm NewA}\theta)\approx qRd\theta/v_{\|}$. The transit average is over a full bounce for trapped ions and over a complete poloidal circuit for the passing ones. The sign of $v_{\|}$ changes at turning points, while $d\tau>0$, giving $\overline{v_{\|}}=0$ for the trapped. For the passing $\overline{v_{\|}}\rightarrow\langle 1/v_{\|}\rangle^{-1}$ in the large aspect ratio limit, where

$\langle X\rangle=[\oint Xd\theta/(\mathbf{B}_{0}\cdot{\rm NewA}\theta)][\oint d\theta/(\mathbf{B}_{0}\cdot{\rm NewA}\theta)]^{-1}$ denotes a flux surface average. It is easy to show that $\langle\int d^{3}v\,f_{0}X\rangle=\langle\int d^{3}v\,f_{0}\overline{X}\rangle$. More details are given in Appendix A.

Equation (29) generalizes the usual electrostatic result to include electromagnetic effects through $A_{\|}$ and $B_{\|}$. It is important to realize that $\Phi$, $A_{\|}$ and $B_{\|}$ are allowed to be slow functions of $\theta$ and $\psi$. Solving (29) by integrating from $t=0$ gives

[TABLE]

where $\tilde{\Phi}^{(0)}=\tilde{\Phi}(t=0)$, $\tilde{A}_{\|}^{(0)}=\tilde{A}_{\|}(t=0)$ and $\tilde{B}_{\|}^{(0)}=\tilde{B}_{\|}(t=0)$ are allowed to be flux functions, $J_{0}\equiv J_{0}(z)$ and $J_{1}\equiv J_{1}(z)$, and $z=k_{\perp}v_{\perp}/\Omega$. The strong poloidal variation in (30) is due to $v_{\|}$ and $Q\propto v_{\|}$, while the poloidal variation of $B_{0}$ is weaker and sometimes unimportant. The electron response $\overline{\tilde{h}_{e}}$ is given by an equation similar to (30), but with $Z\rightarrow-1$, $M\rightarrow m$ and $f_{0}\rightarrow f_{0e}$. We will omit species subscripts to streamline notation except where they are needed to avoid confusion.

Next we form the perturbed quasineutrality equation

[TABLE]

and the two components of Ampère’s law

[TABLE]

and

[TABLE]

where the integrals are taken at fixed $\psi$, $\sum$ denotes a sum over ions and electrons. Notice that $\mathbf{J}_{1}\cdot\mathbf{e}_{\perp}\simeq 0$ since ${\rm NewA}\cdot\mathbf{J}_{1}=0={\rm NewA}\cdot\mathbf{A}$. Using $f_{1}=\tilde{h}_{\ast}\exp[iS(\psi_{\ast})]-(Zef_{0}/T)\tilde{\Phi}\exp[iS(\psi)]$ and $\int d^{3}vh\left\{1,v_{\|},v_{\perp}\sin\varphi\right\}=\int d^{3}v\tilde{h}_{\ast}\left\{1,v_{\|},v_{\perp}\sin\varphi\right\}\exp[i(L-Q+S)]=\int d^{3}v\tilde{h}_{\ast}\left\{J_{0},v_{\|}J_{0},iv_{\perp}J_{1}\right\}\exp[i(S-Q)]$, equations (31), (32) and (33) give

[TABLE]

[TABLE]

and

[TABLE]

Note that the perpendicular Ampère’s law (33) has become perpendicular pressure balance (36).

Equation (30) represents the formal solution of an initial value problem for a time scale long compared to the periodic gyromotion about guiding centers. It depends on the initial value $\tilde{h}_{\ast}(0)$, which we are free to choose arbitrarily as long as it satisfies Maxwell’s equations. Our particular choice for $\tilde{h}_{\ast}(0)$ will be motivated in the next section.

4 Choice of initial condition

We must pick $\tilde{h}_{\ast}$ to satisfy the Maxwell equations at $t=0$, while also obtaining a convenient and sensible form for $\overline{\tilde{h}_{\ast}}$ to evaluate the long time relaxation behavior of the zonal flow response. The non-transit averaged initial condition for $\tilde{h}_{\ast}(t=0)$ must depend on $\theta$ as well as $\psi_{\ast}$, $v$, $\mu$ to generate a transiently evolving zonal flow (the GAM) or else it will be a homogeneous solution to the non-transit averaged Vlasov equation. However, we do not allow $\tilde{h}_{\ast}(t=0)$ to depend on gyrophase since we assume any initial transient response associated with gyromotion has already damped to its classical polarization level.

We require Maxwell’s equations to be satisfied at $t=0$, which we take to mean after many gyrations, but much less than the time for a poloidal bounce or transit to be completed. Consequently, at $t=0$ we must satisfy (34)-(36).

A GAM develops on a time scale of the order of a transit or bounce time and much longer than a gyration period. It oscillates as it damps away to the residual zonal flow level. The initial conditions are sometimes viewed as approximating the turbulent sources (Rosenbluth & Hinton, 1998; Sugama & Watanabe, 2005) of charge and current densities on a time much less than a transit or bounce time, but after many gyrations. When the long time behavior of the system is studied electrostatically we may use the transit average result

[TABLE]

where $\overline{\tilde{h}_{\ast}(t)}$ is given by the electrostatic limit of (30).

Next we explain how we choose $\tilde{h}_{\ast}(0)$ to include magnetic perturbations. For our results to have the required generality we do not flux surface average quasineutrality and the two components of Ampère’s law. Moreover, we desire forms at $t=0$ in which $\tilde{\Phi}^{(0)}$, $\tilde{A}_{\|}^{(0)}$ and $\tilde{B}_{\|}^{(0)}$ terms only contribute to quasineutrality, and the parallel and perpendicular components of Ampère’s law, respectively. To avoid the need for complicated velocity space structure, we do not consider arbitrary wavenumbers (Xiao et al., 2007a; Xiao & Catto, 2006b). We do manipulations consistent with an expansion in $k_{\perp}v_{\perp}/\Omega\ll 1$, but to avoid lengthy expressions, we prefer to keep the finite Larmor radius terms in the form of Bessel functions for a while and make the expansions explicit later. To simplify our treatment and properly recover the $k_{\perp}v_{\perp}/\Omega\ll Q\ll 1$ limit electrostatically we can use the simple form

[TABLE]

where since we are only interested in $k_{\perp}v_{\perp}/\Omega\ll 1$ the zeroes of the Bessel function are of no concern. A nice discussion of the difference between treating electrostatic zonal flows as an initial value problem rather than as a turbulent source in quasineutrality is presented in Sec. 3 of Monreal et al. (2016). Moreover, their paper and the references therein should be consulted to understand the differences between the residual zonal flow behaviour in stellarators and tokamaks.

The preceding expression motivates us to assume

[TABLE]

where we have introduced the species dependent quantities

[TABLE]

and

[TABLE]

with $z=k_{\perp}v_{\perp}/\Omega$, $2J_{1}/zJ_{0}\approx 1+z^{2}/8$ giving $\alpha\approx k_{\perp}^{2}T/2M\Omega^{2}\ll 1$ for $z\ll 1$, and the species independent quantity

[TABLE]

In the preceding $\beta=8\pi\sum nT/B_{0}^{2}=8\pi(n_{i}T_{i}+n_{e}T_{e})/B_{0}^{2}=\beta_{i}+\beta_{e}$ and $\omega_{p}^{2}=4\pi\sum Z^{2}e^{2}n/M=4\pi e^{2}[(Z^{2}n_{i}/M)+(n_{e}/m)]=\omega_{pi}^{2}+\omega_{pe}^{2}$. We will also make use of the definitions $v_{i}=(2T_{i}/M)^{1/2}$, $\rho_{i}=v_{i}/\Omega_{i}$, $\Omega_{i}=ZeB_{0}/(Mc)$, and $\rho_{pi}=q\rho_{i}/\epsilon$.

The functional form of $H$ in (39) is chosen so at $t=0$ it does not alter quasineutrality (34) since

[TABLE]

Moreover the $\alpha$ in $H$ is chosen so $\tilde{\Phi}^{(0)}$ will not enter perpendicular Ampère’s law by taking

[TABLE]

The $\tilde{A}_{\|}^{(0)}$ term does not enter quasineutrality because its integral is odd in $v_{\|}$, and $\tilde{B}_{\|}^{(0)}$ does not enter because $u$ and $\beta$ are species independent as we can use unperturbed quasineutrality. We can also see that parallel Ampère’s law, (35), is satisfied at $t=0$ since the $\tilde{\Phi}^{(0)}$ and $\tilde{B}_{\|}^{(0)}$ integrals are odd in $v_{\|}$ leaving

[TABLE]

To satisfy the perpendicular Ampère’s law, (36), at $t=0$ we require (40) to (42) to be satisfied as can be seen from

[TABLE]

where we use unperturbed quasineutrality as well as $\int d^{3}vf_{0}v_{\perp}^{2}[(2J_{1}H/zJ_{0})-1]=0$.

Apart from the need to satisfy Maxwell’s equations, $\tilde{h}_{\ast}(t=0)$ is arbitrary. Other choices for $\tilde{h}_{\ast}(t=0)$ may be used, but (39) is sufficient for our purpose. It has the transit average

[TABLE]

where the weak $\theta$ dependences of $u$ and $\alpha$ are neglected.

Before closing this section we prove that $\tilde{\Phi}$ is a lowest order flux function for fully electromagnetic initial conditions. For this demonstration we will only be interested in the $L\sim k_{\perp}v_{\perp}/\Omega\ll Q\sim k_{\perp}v_{\|}/\Omega_{p}\ll 1$ limit. Consequently, to lowest order we let $k_{\perp}\rightarrow 0$ to make $e^{iQ}\rightarrow 1$, $J_{0}\rightarrow 1$, $J_{1}\rightarrow k_{\perp}v_{\perp}/2\Omega$, $\alpha\rightarrow 0$, $H\rightarrow 1$ and $u\rightarrow 1$ in (34)-(36), and then use $\tilde{h}_{\ast}\approx\overline{\tilde{h}_{\ast}}$ to obtain quasineutrality in the form

[TABLE]

with

[TABLE]

and

[TABLE]

where we have set $\alpha=0$ and $u=1$ since we let $k_{\perp}\rightarrow 0$ and used $\overline{v_{\perp}^{2}/2B_{0}}=\mu=v_{\perp}^{2}/2B_{0}$. The integrals odd in $v_{\|}$ do not contribute to quasineutrality, of course. In addition, unperturbed quasineutrality prevents the $\tilde{B}_{\|}^{(0)}$ and $\tilde{B}_{\|}$ from contributing to perturbed quasineutrality as well since they result in the integrals $\sum Ze\int d^{3}vf_{0}\overline{B_{0}}=(B_{0}/2\underline{B}_{0})(\int_{0}^{\underline{B}_{0}/B_{0}}d\lambda\overline{B_{0}}/\xi)\sum Zen=0$ and

[TABLE]

where $B_{0}/\underline{B}_{0}\simeq R_{0}(\psi)/R=1-\epsilon(\psi)\cos\theta+\mathcal{O}(\epsilon^{2})$, $\underline{B}_{0}=\langle B_{0}\rangle$, $\epsilon=r(\psi)/R_{0}(\psi)$, $\lambda=2\mu\underline{B}_{0}/v^{2}$, $\xi=\sqrt{1-\lambda B_{0}/\underline{B}_{0}}$, and $d^{3}v\rightarrow\sum_{\mathrm{sgn}v_{\|}}\pi B_{0}v^{3}dvd\lambda/\underline{B}_{0}|v_{\|}|$ (where $\sum_{\mathrm{sgn}\,v_{\|}}$ can be replaced by a factor 2 when the integrand is even in $\mathrm{sgn}\,v_{\|}$). Consequently, assuming $k_{\perp}v_{\perp}/\Omega\ll Q\ll 1$ the lowest order perturbed quasineutrality equation becomes

[TABLE]

which requires

[TABLE]

and is satisfied if $\tilde{\Phi}=\langle\tilde{\Phi}\rangle$. Therefore, we may safely assume $\tilde{\Phi}$ is a flux function to lowest order. To see this more rigorously, we multiply (52) by $\tilde{\Phi}$ to form $\tilde{\Phi}^{2}\int d^{3}vf_{0}=\int d^{3}vf_{0}\tilde{\Phi}^{2}=\tilde{\Phi}\int d^{3}vf_{0}\overline{\tilde{\Phi}}=\int d^{3}vf_{0}\overline{\tilde{\Phi}}\tilde{\Phi}$, then flux surface average to obtain

[TABLE]

As a result, using $\left\langle\int d^{3}vf_{0}\tilde{\Phi}^{2}\right\rangle=\left\langle\int d^{3}vf_{0}\overline{\tilde{\Phi}^{2}}\right\rangle$ and $\left\langle\int d^{3}vf_{0}\overline{\tilde{\Phi}}\tilde{\Phi}\right\rangle=\left\langle\int d^{3}vf_{0}\overline{\tilde{\Phi}}^{2}\right\rangle$ we find

[TABLE]

Therefore, we need

[TABLE]

but $\tilde{\Phi}=\tilde{\Phi}(\psi,\theta)$, while $\overline{\tilde{\Phi}}=\overline{\tilde{\Phi}}(\psi,\lambda)$, so it must be that

[TABLE]

to lowest order.

To generalize our results for quasineutrality and to treat both components of Ampère’s law, we extend our initial condition (50) to include finite orbit effects by employing (47) in (30).

In the calculation of the electrostatic zonal flow residual (Rosenbluth & Hinton, 1998) the non-adiabatic electron response could be neglected as small in $(m/M)^{1/2}$, where $m$ is the electron mass, however here it sometimes needs to be retained, since the electromagnetic terms in the non-adiabatic response are proportional to the thermal speed of the species. In the next section we will use the preceding results to form quasineutrality and Ampère’s law.

5 Quasineutrality and Ampère’s law

In this section we form and consider the non-flux surfaced average components of Ampère’s law to demonstrate that poloidal variation of the parallel vector potential and the parallel magnetic field must be retained. Indeed, we will find that there are cases for which these field responses have strong poloidal variation. In addition, we will perform a more complete evaluation of quasineutrality once we have examined the two components of Ampère’s law.

To perform the derivation of the two components of Ampère’s law we must realize that $\tilde{A}_{\|}$ and $\tilde{B}_{\|}$ are not normally flux functions. Indeed, even when $\tilde{B}_{||}^{(0)}=0=\tilde{A}_{||}^{(0)}$ we will find they both have poloidal variation. Retaining $e^{iQ}$ modifications, but ignoring $L\sim k_{\perp}v_{\perp}/\Omega$ corrections as unimportant except in the $\tilde{\Phi}$ and ${\tilde{\Phi}^{(0)}}$ terms, we use

[TABLE]

with

[TABLE]

5.1 Ampère’s law

Inserting the preceding into

[TABLE]

and

[TABLE]

we obtain

[TABLE]

and

[TABLE]

Recalling that $\tilde{A}_{\|}$ and ${\tilde{B}_{||}}$ need not be flux functions and that the unperturbed magnetic field gives rise to $\cos\theta$ dependence, we Fourier decompose and retain only the leading poloidal dependence by writing

[TABLE]

and

[TABLE]

with the coefficients $\langle\tilde{A}_{\|}\rangle$, $\langle{\tilde{B}_{||}}\rangle$, $\langle\tilde{A}\rangle$, $\langle\tilde{B}\rangle$, $a$, and $b$ flux functions, and $\mathbf{B}_{0}\cdot{\rm NewA}\theta\propto 1-(\epsilon-\Delta_{S}^{\prime})\cos\theta$ in the flux surface averages for our Shafranov shifted circular flux surface model, where $\Delta_{S}^{\prime}=d\Delta_{S}/dr\sim\epsilon$ . We do not assume an ordering of $a$ and $b$ relative to $\langle\tilde{A}_{\|}\rangle$ and $\langle{\tilde{B}_{||}}\rangle$ as we will calculate all four of these coefficients. In addition, we recall the electrostatic result of Rosenbluth & Hinton (1998),

[TABLE]

with $\gamma\approx 1.64$. The preceding allows us to assume ${q^{2}}{\epsilon^{-1/2}}\left\langle\tilde{\Phi}\right\rangle\sim{\tilde{\Phi}^{(0)}}$ for $\beta\ll 1$ and $\epsilon\ll 1$.

Expanding for $Q=k_{\perp}v_{\|}/\Omega_{p}\ll 1$, inserting $H$, neglecting $\epsilon$ corrections to $k_{\perp}^{2}\rho_{i}^{2}$ terms, ignoring all $Q^{2}$, $\beta$ and $\epsilon^{3/2}$ corrections to the ${\tilde{B}_{||}}$ and $\tilde{B}_{||}^{(0)}$ terms in (62), recalling $\tilde{\Phi}=\left\langle\tilde{\Phi}\right\rangle$, and using unperturbed quasineutrality, perpendicular Ampère’s law becomes

[TABLE]

Here and hereafter, we use $Q\propto v_{\|}/B_{0}$ and then treat $k_{\perp}$ as a flux function, except in $k_{\perp}^{2}\tilde{A}_{||}^{(0)}$ terms from (59) and $k_{\perp}^{2}\tilde{A}_{||}$ from (61).

In (63) we neglect $Q$ corrections for $\tilde{A}_{\|}$ and $\tilde{A}_{||}^{(0)}$ terms, insert $H$, and continue to use $\tilde{\Phi}=\left\langle\tilde{\Phi}\right\rangle$, to find that parallel Ampère’s law is

[TABLE]

where we used (140), (148), and (155) and similar integrals for $\left\langle\tilde{\Phi}\right\rangle$ and ${\tilde{\Phi}^{(0)}}$ terms to see that due to unperturbed quasineutrality$\sum({Z^{2}}{e^{2}}/T)\int{d^{3}}v{f_{0}}{v_{||}}(Q-\overline{Q})=0$. The zonal flow responses of the perturbed parallel magnetic field and vector potential in (67) and (68) are generated by the polarizations terms associated with poloidally varying departures from flux surfaces.

Next, we flux surface average using $\left\langle\int{d^{3}}v{f_{0}}\overline{X}Y\right\rangle=\left\langle\int{d^{3}}v{f_{0}}X\overline{Y}\right\rangle$, to obtain the flux surface averaged perpendicular Ampère’s law

[TABLE]

where we used (162), (185), (188), (189), and (194). The $B_{0}$ factors in (67) must be treated with care, even though $\overline{v_{\|}\cos\theta}=0$ for the trapped and $\overline{v_{\|}\cos\theta}\sim\epsilon v$ for the passing, since they alter the coefficient. For now we retain the $\epsilon^{1/2}$ smaller $\left\langle\tilde{\Phi}\right\rangle$ and ${\tilde{\Phi}^{(0)}}$ terms compared to the $\gamma{\epsilon^{-1/2}}{q^{2}}\left\langle\tilde{\Phi}\right\rangle$ term.

In the parallel Ampère’s law for now we keep $\mathcal{O}(\epsilon^{3/2})$ corrections to the typically large $\omega_{p}^{2}/k_{\perp}^{2}{c^{2}}$ terms obtained by using (148). By keeping the $\epsilon^{3/2}$ correction to $\tilde{A}_{\|}$ terms we are making the easily satisfied assumption that ${\epsilon^{3/2}}\omega_{p}^{2}/k_{\perp}^{2}{c^{2}}\gg{q^{2}}\beta$, or roughly $M/m\gg\epsilon^{1/2}k_{\perp}^{2}\rho_{pi}^{2}$, since we already neglected the more complicated $Q^{2}$ corrections to the same terms. The flux surface averaged parallel Ampère’s law then becomes

[TABLE]

where the constant $\sigma\approx 5.3$ is evaluated in Appendix A, and we make use of (148), (153), (155), (162), (195) and (202). We use (64) and keep the Shafranov shift, $\Delta_{S}$, terms from $k_{\perp}\propto RB_{p}\propto 1-\Delta_{S}^{\prime}\cos\theta$ in the $\langle\tilde{A}_{\|}\rangle$ and $\tilde{A}_{\|}^{(0)}$ terms that arise from (59) and the left hand side of (61). To simplify the flux surface averaged parallel Ampère’s law we ignored $\beta{\tilde{B}_{||}}$ and $\beta\tilde{B}_{||}^{(0)}$ terms. Also, we have neglected all $b$ terms because they are multiplied by $\beta$. Based on (66) we must retain the $\sigma{q^{2}}{\epsilon^{-1/2}}\left\langle\tilde{\Phi}\right\rangle$ term in (70), but will ignore the $\epsilon$ and ${\epsilon^{3/2}}$ smaller term $(5\gamma{\epsilon^{1/2}}+4\epsilon)\left\langle\tilde{\Phi}\right\rangle$, as well as the $(\gamma{\epsilon^{1/2}}+\epsilon){\tilde{\Phi}^{(0)}}$ term, all associated with the Bessel function corrections.

The terms $\gamma$ and $\sigma$ in (69) and (70) arise from the trapped particle responses to any initial perturbation, while the $a\epsilon$ and $a\Delta_{S}^{\prime}$ terms are passing responses. We estimate the different behaviour of the passing and trapped (and barely passing) by using $\overline{{v_{||}}}=0$ and estimating ${v_{||}}\sim{\epsilon^{1/2}}v$ and ${\int d^{3}}v\propto{\epsilon^{1/2}}$ for the trapped particles, while using $\overline{{v_{||}}}\simeq{v_{||}}\left[1+\mathcal{O}(\epsilon)\right]$ for the passing ones.

Next, we subtract the flux surface averaged equations from the full equations. For the perpendicular Ampère’s law we use (185), (188), (189), and (194) to find

[TABLE]

where Bessel function terms are ignored since any $\cos\theta$ dependence is smaller by $\epsilon$, and the constants $\chi\approx 0.11$ and $\varsigma\approx 5.3$ are evaluated in Appendix A. The $\left\langle\tilde{A}_{\|}\right\rangle$, $\tilde{A}_{||}^{(0)}$, and $\left\langle\tilde{\Phi}\right\rangle$ terms are from passing contributions, while the ${\epsilon^{3/2}}a$ term is a trapped contribution.

To form the difference equation for the parallel Ampère’s law we ignore $\beta\left\langle{\tilde{B}_{||}}\right\rangle$, $\beta\tilde{B}_{||}^{(0)}$ and $\beta b$ terms. As a result, the poloidally varying parallel Ampère’s law reduces to

[TABLE]

where we make use of use (148), (153), (155), (162), (195), and (202), and account for $k_{\perp}\propto 1-\Delta_{S}^{\prime}\cos\theta$. Only passing contributions enter for $\left\langle\tilde{A}_{\|}\right\rangle$, $\tilde{A}_{||}^{(0)}$, $a$, and $\tilde{B}_{||}^{(0)}$ terms, and only Bessel terms enter for ${\tilde{\Phi}^{(0)}}$, while for $\left\langle\tilde{\Phi}\right\rangle$ both trapped and Bessel contributions enter.

5.2 Parallel Ampère’s law limit

Simulations are sometimes run ignoring the perpendicular Ampère’s law. If we do the same by dropping all $b$, and $\tilde{B}_{||}^{(0)}$ terms in parallel Ampère’s law, then (72) becomes

[TABLE]

and (70) reduces to

[TABLE]

where we have used $\beta\ll 1$ so we can use (66) to eliminate $\left\langle\tilde{\Phi}\right\rangle$.

If $\beta\rightarrow 0$ at finite $k_{\perp}$ so that $\omega_{p}^{2}/k_{\perp}^{2}{c^{2}}\to 0$ then (73) and (74) reduce to the large skin depth results

[TABLE]

and

[TABLE]

where we neglect ${\epsilon^{1/2}}$ corrections to ${\tilde{\Phi}^{(0)}}$ terms. Interestingly, for $\tilde{A}_{||}^{(0)}=0$ we see that $\left\langle\tilde{A}_{\|}\right\rangle/a\simeq\sigma/2\gamma\simeq\mathcal{O}(1)$, meaning that substantial poloidal variation occurs when ${\tilde{\Phi}^{(0)}}\neq 0$. When ${\tilde{\Phi}^{(0)}}=0$, $\left\langle\tilde{A}_{\|}\right\rangle\simeq\tilde{A}_{||}^{(0)}(1-\gamma{\epsilon^{3/2}})$ and $a\simeq-\epsilon\left\langle\tilde{A}_{\|}\right\rangle$, and poloidal variation is somewhat weak.

More interestingly, we consider finite $\beta\ll 1$ by allowing $\epsilon^{3/2}\gg k_{\perp}^{2}{c^{2}}/\omega_{p}^{2}$ or ${q^{2}}{\beta_{i}}\epsilon^{1/2}\gg Zmk_{\perp}^{2}\rho_{pi}^{2}/M$. Continuing to neglect ${\epsilon^{1/2}}$ corrections in this small skin depth limit we find

[TABLE]

and

[TABLE]

When ${\tilde{\Phi}^{(0)}}=0$, $\left\langle\tilde{A}_{\|}\right\rangle\simeq\tilde{A}_{||}^{(0)}\simeq-a/\epsilon$ so only weak poloidal variation occurs. However, for $\tilde{A}_{||}^{(0)}\simeq 0$, we find that $a/\left\langle\tilde{A}_{\|}\right\rangle\simeq 2$, which results in strong poloidal variation.

In both the small and large skin depth limits, when ${\tilde{\Phi}^{(0)}}=0$, only small changes in $\tilde{A}_{\|}$ occur since $\tilde{A}_{\|}\simeq\left\langle\tilde{A}_{\|}\right\rangle\simeq\tilde{A}_{||}^{(0)}\gg a$. However, for $\tilde{A}_{||}^{(0)}\simeq 0$, we see that poloidal variation in $\tilde{A}_{\|}$ arises due to ${\tilde{\Phi}^{(0)}}$ with $a/\left\langle\tilde{A}_{\|}\right\rangle\sim 1$ in both limits.

5.3 Full Ampère’s law

The general case retains the perpendicular Ampère’s law. Using (70) and (72) with (66) inserted and ignoring ${\epsilon^{1/2}}$ corrections we find

[TABLE]

and

[TABLE]

In addition, (69) and (71) reduce to

[TABLE]

and

[TABLE]

where we used the lowest order version of (66),

[TABLE]

When $\tilde{B}_{||}^{(0)}=0$ we may employ (75)-(78) to find $\left\langle{\tilde{B}_{||}}\right\rangle$ and $b$ from (81) and (82). In the $\omega_{p}^{2}/k_{\perp}^{2}{c^{2}}\to 0$ limit (75) and (76) are used along with (81) and (82) to find

[TABLE]

and

[TABLE]

Consequently, when ${\tilde{\Phi}^{(0)}}=0$ we see that $\left\langle{\tilde{B}_{||}}\right\rangle/b\simeq-5\gamma{\epsilon^{1/2}}/6$ so strong poloidal variation occurs in ${\tilde{B}_{||}}$. For $\tilde{A}_{||}^{(0)}=0$ the poloidal variation of ${\tilde{B}_{||}}$ is strong since $b/\left\langle{\tilde{B}_{||}}\right\rangle\simeq{\varsigma}/\gamma\approx 3.2$.

For $\tilde{B}_{\|}^{(0)}=0$, $\beta\ll 1$, and $\epsilon^{3/2}\gg k_{\perp}^{2}{c^{2}}/\omega_{p}^{2}$ we use (77) and (78) to find the small skin depth forms

[TABLE]

and

[TABLE]

When $\tilde{A}_{||}^{(0)}=0$ the poloidal variation of ${\tilde{B}_{||}}$ is strong since $b/\left\langle{\tilde{B}_{||}}\right\rangle\simeq{\varsigma}/\gamma\approx 3.2$. For ${\tilde{\Phi}^{(0)}}=0$ it is also strong with $b/\left\langle{\tilde{B}_{||}}\right\rangle\simeq-[1+(2k_{\perp}^{2}c^{2}/\epsilon^{2}\omega_{p}^{2})]$. Consequently, in the small skin depth limit, for $\tilde{B}_{||}^{(0)}=0$ the poloidal variation of ${\tilde{B}_{||}}$ is strong when $\tilde{A}_{||}^{(0)}=0$ and weak when ${\tilde{\Phi}^{(0)}}=0$.

Next we consider the case ${\tilde{\Phi}^{(0)}}=0$ and $\tilde{A}_{||}^{(0)}=0$ when $\tilde{B}_{||}^{(0)}\neq 0$ by using (79)-(82). The $\omega_{p}^{2}/k_{\perp}^{2}{c^{2}}\to 0$ limit gives

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Consequently, the poloidal variation of ${\tilde{B}_{||}}$ is weak (with $b/\left\langle{\tilde{B}_{||}}\right\rangle\simeq\epsilon$), while the poloidal variation of $\tilde{A}_{\|}$ is strong (with $\left\langle\tilde{A}_{\|}\right\rangle/a\simeq\gamma{\epsilon^{1/2}}$) in this large skin depth limit.

In the small skin depth limit $\epsilon^{3/2}\gg k_{\perp}^{2}{c^{2}}/\omega_{p}^{2}$ for finite $\beta\ll 1$ and $\tilde{\Phi}^{(0)}=0=\tilde{A}_{\|}^{(0)}$ we use (79) and (80) to find $\epsilon a\simeq 2(1+\gamma\epsilon^{3/2})\langle\tilde{A}_{\|}\rangle$. Then (79) to (82) give

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

As a result, for this $\tilde{B}_{||}^{(0)}\neq 0$ case with ${\tilde{\Phi}^{(0)}}=0$ and $\tilde{A}_{||}^{(0)}=0$, the poloidal variation of $\tilde{A}_{\|}$ is strong with $\left\langle\tilde{A}_{\|}\right\rangle/a\simeq\epsilon/2$, while the poloidal variation of ${\tilde{B}_{||}}$ varies from weak to strong with $b/\left\langle{\tilde{B}_{||}}\right\rangle\simeq\epsilon$ to $1$.

5.4 Quasineutrality

To complete our treatment of the Maxwell equations we need to form quasineutrality with its finite $\beta$ effects retained using

[TABLE]

After inserting (30) and (47), we only retain Bessel function modifications to $\tilde{\Phi}\simeq\left\langle\tilde{\Phi}\right\rangle$ and ${\tilde{\Phi}^{(0)}}$ terms. We then find

[TABLE]

Expanding the Bessel functions and the exponentials in $Q$, inserting $H$ with $\alpha=k_{\perp}^{2}T/2M{\Omega^{2}}$, and using unperturbed quasineutrality leaves

[TABLE]

We have already shown that the poloidal variation of $\tilde{\Phi}$ is negligible so we only require the flux surface average of (98), which is

[TABLE]

where we neglect $\beta\left\langle{\tilde{B}_{||}}\right\rangle$, $\beta\tilde{B}_{||}^{(0)}$ and $\beta b$ corrections.

Performing the integrals using (148), (198), (201), (204), and (205), and noting that only the trapped ion and Bessel contributions matter we obtain

[TABLE]

In the limit in which we take $\beta\rightarrow 0$ at finite $k_{\perp}$ (such that $\omega_{p}^{2}/k_{\perp}^{2}{c^{2}}\to 0$), (100) with (76) and (77) inserted gives the large skin depth expression

[TABLE]

In the more interesting finite limit, for which $\epsilon^{3/2}\gg k_{\perp}^{2}{c^{2}}/\omega_{p}^{2}$, using (77), (78), (92), and (93) we find the small skin depth result

[TABLE]

Clearly, the modification of the electrostatic limit is small when $\tilde{B}_{||}^{(0)}=0=\tilde{A}_{||}^{(0)}$. From (102), we see that initial magnetic perturbations as large as ${\epsilon^{1/2}}q{\rho_{pi}}{v_{i}}\tilde{B}_{||}^{(0)}/c\sim\beta{\tilde{\Phi}^{(0)}}$ and/or ${\epsilon^{3/2}}{q^{2}}{k_{\perp}}{\rho_{pi}}{v_{i}}\tilde{A}_{||}^{(0)}/c\sim{\tilde{\Phi}^{(0)}}$ are required to obtain order unity corrections to the electrostatic response. For the same size initial perturbations in the expression (101) the response to $\tilde{B}_{||}^{(0)}$ remains the same, but that due to $\tilde{A}_{||}^{(0)}$ is very small.

6 Zonal flow responses in terms of initial field values: 15 test cases

In this section we present a summary of our results in the small skin depth limit. Before doing so we note that

[TABLE]

with $\beta={\beta_{i}}(1+\tau)$ and $\tau=Z{T_{e}}/{T_{i}}$. We, of course, neglect mass ratio corrections in the plasma frequency.

The large skin depth limit of $\beta\rightarrow 0$ requires ${q^{2}}\beta\ll Zmk_{\perp}^{2}\rho_{pi}^{2}/M\ll Zm/M$ since we must keep $k_{\perp}^{2}\rho_{pi}^{2}\ll 1$ for our analysis to hold. Consequently, the more interesting limit is that of small skin depth ${\epsilon^{2}}\sim\epsilon{\Delta^{\prime}_{S}}\gg k_{\perp}^{2}{c^{2}}/\omega_{p}^{2}$ and small ${\beta_{i}}$ ($1\gg\beta_{i}\gtrsim Zmk_{\perp}^{2}\rho_{pi}^{2}/q^{2}M$), which will be our focus. In the following we summarize our results for this limit. We will find that electromagnetic perturbations give an electrostatic potential response comparable to the electrostatic limit of $\left\langle\tilde{\Phi}\right\rangle\simeq{\epsilon^{1/2}}{\tilde{\Phi}^{(0)}}/\gamma{q^{2}}$ when ${v_{i}}\tilde{A}_{\|}^{(0)}{\epsilon^{3/2}}{q^{2}}{k_{\perp}}{\rho_{pi}}/c\sim{\tilde{\Phi}^{(0)}}$ or ${\epsilon^{1/2}}q{\rho_{pi}}{v_{i}}\tilde{B}_{||}^{(0)}/c\sim\beta{\tilde{\Phi}^{(0)}}$.

6.1 Electromagnetic response to an electrostatic initial perturbation ($\tilde{B}_{||}^{(0)}=0=\tilde{A}_{\|}^{(0)}$)

In this case $\left\langle\tilde{\Phi}\right\rangle/{\tilde{\Phi}^{(0)}}$ has a very small linear in $\beta$ correction to the electrostatic Rosenbluth and Hinton form as seen from (102),

[TABLE]

where $\sigma\approx 5.3$. Moreover, $\left\langle\tilde{B}_{\|}\right\rangle{\rho_{i}}{v_{i}}/c{\tilde{\Phi}^{(0)}}$ and $b{\rho_{i}}{v_{i}}/c{\tilde{\Phi}^{(0)}}$ are linear in $\beta$ and are seen from (86) and (87) to be given by

[TABLE]

and

[TABLE]

From (77) and (78) we find

[TABLE]

and

[TABLE]

where all electromagnetic responses are proportional to $\beta$. The correction to the electrostatic result is very small, but all electromagnetic responses are larger and of the same order so strong poloidal variation occurs.

6.2 Electromagnetic response to a plucked field line initial condition ($\tilde{B}_{||}^{(0)}=0={\tilde{\Phi}^{(0)}}$)

Using (102) gives the electrostatic potential response for this $\tilde{A}_{\|}^{(0)}\neq 0$ case to be

[TABLE]

For this case, the parallel vector potential responses from (77) and (78) are given by

[TABLE]

and

[TABLE]

showing that $\tilde{A}_{\|}$ will only be very slightly perturbed from $\tilde{A}_{\|}^{(0)}$ with weak poloidal variation due to field line plucking satisfying $R\tilde{A}_{\|}\simeq R_{0}(\psi)\tilde{A}_{\|}^{(0)}$. The $\beta$ proportional $\tilde{B}_{\|}$ responses from (86) and (87) give $b$ and $\left\langle{{\tilde{B}}_{||}}\right\rangle$ as comparable as seen from

[TABLE]

and

[TABLE]

6.3 Electromagnetic response to a compressed field line

initial condition ($\tilde{A}_{\|}^{(0)}=0={\tilde{\Phi}^{(0)}}$)

The response $c{k_{\perp}}\left\langle\tilde{\Phi}\right\rangle/{v_{i}}\tilde{B}_{||}^{(0)}$ is proportional to $1/\beta$ as seen from (102) so in this $\tilde{B}_{||}^{(0)}\neq 0$ case we multiply through by $\beta$ to form

[TABLE]

The responses $\left\langle\tilde{B}_{\|}\right\rangle/\tilde{B}_{||}^{(0)}$ and $b/\tilde{B}_{||}^{(0)}$ are given by (94) and (95) for this field line stretching case:

[TABLE]

and

[TABLE]

The responses ${k_{\perp}}\left\langle\tilde{A}_{\|}\right\rangle/i\tilde{B}_{||}^{(0)}$ and ${k_{\perp}}a/i\tilde{B}_{||}^{(0)}$ are obtained from (92) and (93)

[TABLE]

and

[TABLE]

Stretching or compressing field lines causes very strong poloidal variation in $\tilde{A}_{\|}$, but weak poloidal variation in $\tilde{B}_{\|}$ satisfying $R_{0}(\psi)\tilde{B}_{\|}\simeq R\tilde{B}_{\|}^{(0)}$.

7 Conclusions

We have derived approximate analytical expressions for long wavelength, collisionless zonal flow residual responses at low $\beta$ in the Shafranov shifted, circular cross section, large aspect ratio limit of a tokamak. To do so, we formulate and solve a Maxwell-Vlasov description in the form of an initial value problem, retaining the fully self-consistent spatial perturbations in the electric and magnetic fields. As zonal flow perturbations are axisymmetric, the only magnetic field topology change allowed is the switch between rational and irrational field lines within a flux surface - no magnetic islands can be formed. The choice of the initial condition in the non-adiabatic part of the distribution function must be consistent with Maxwell’s equations, but it is otherwise arbitrary. The specific choice we make is motivated by the desire to recover the usual long wavelength result in the electrostatic limit that has the residual proportional to the ratio of the classical polarization over the classical plus neoclassical polarization. Also, our choice of initial conditions is such that at t = 0 the initial electrostatic, and shear- and compressional magnetic perturbations contribute only to quasineutrality, the parallel and the perpendicular components of Ampère’s law, respectively. This form is convenient since then the initial conditions in the amplitude of the fields can be chosen independently.

The results we obtain are expected to prove useful for testing turbulent electromagnetic gyrokinetic codes just as Rosenbluth and Hinton has proven useful in the electrostatic limit. The electromagnetic case is of course far more complex, and further complicated by the fact that the parallel vector potential and the parallel magnetic field responses are no longer simply flux functions. Our electromagnetic zonal flow responses provide 15 meaningful tests of fully electromagnetic gyrokinetic turbulence codes. We focus on the small skin depth limit in this section and sections 5 and 6, but section 5 also gives large skin depth results.

For a pure electrostatic initial condition (${\tilde{\Phi}^{(0)}}\neq 0$, $\tilde{B}_{||}^{(0)}=0=\tilde{A}_{\|}^{(0)}$), the usual long wavelength Rosenbluth and Hinton electrostatic zonal flow result is recovered with only a small $\beta$ correction as seen from (104). However, the responses of the parallel vector potential and parallel magnetic fields will have $\cos\theta$ dependence as well flux surface averaged responses and are given by (105)-(108). Indeed the poloidally varying and flux surface averaged responses are comparable for $\tilde{B}_{\|}$ and $\tilde{A}_{\|}$.

The shear Alfvén field line plucking initial condition ($\tilde{A}_{\|}^{(0)}\neq 0$, $\tilde{B}_{||}^{(0)}=0={\tilde{\Phi}^{(0)}}$), also give stronger flux surface averaged $\tilde{A}_{\|}$ and $\tilde{B}_{\|}$ responses as can be seen from (110)-(113). Because the initial perturbation is electromagnetic, the $\left\langle\tilde{\Phi}\right\rangle$ as well as the $\tilde{A}_{\|}$ responses are free of any $\beta$ multipliers as can be seen from (109)-(111). The compressional Alfvén responses of (112) and (113) contain $\beta$ multipliers. Only weak poloidal variation is found in this case for $\tilde{A}_{\|}$. However, the poloidal variation of $\tilde{B}_{\|}$ is important with $b\sim\langle\tilde{B}_{\|}\rangle$.

When the initial perturbation is pure field compression ($\tilde{B}_{||}^{(0)}\neq 0$, $\tilde{A}_{\|}^{(0)}=0={\tilde{\Phi}^{(0)}}$) the response of $\left\langle\tilde{\Phi}\right\rangle$ appears very large since it is proportional to $1/\beta$ as can be seen from (6.12). In this case $\left\langle\tilde{\Phi}\right\rangle$ is better viewed as an order unity response to a $\tilde{B}_{||}^{(0)}\propto\beta$ initial condition. Then the smallness of $\tilde{B}_{||}^{(0)}$ cancels the $1/\beta$ dependence making $\beta\tilde{B}_{||}^{(0)}$ independent of $\beta$. The $\tilde{B}_{\|}$ response is given by (115) and (116) with poloidal variation again weak. However, the poloidal variation of $\tilde{A}_{\|}$ response is $1/\epsilon$ stronger than the flux surface averaged response as seen from (117) and (118).

**Acknowledgments

**Work supported by U. S. Department of Energy grants at DE-FG02-91ER-54109 at MIT. IP is supported by the Intenational Career Grant (Dnr. 330-2014-631) of Vetenskapsrådet, and Marie Sklodowska Curie Actions, Cofund, Project INCA 600398. This research topic was suggested to the authors by Alex Schekochihin of Oxford in March of 2013 while we all enjoyed the hospitality and support of the Wolfgang Pauli Institute in Vienna, Austria. Fortunately, our friendship has survived in spite of that seemingly harmless suggestion. Indeed, P.J.C. is indebted to Alex for accommodations and hospitality at Merton College during subsequent collaborative visits to the Rudolf Peierls Centre for Theoretical Physics at Oxford University in 2016 and 2017. In addition, the authors collaborated during the Gyrokinetic Theory Working Group Meeting 2016 in Madrid in September.

Appendix A Endless Integrals

For any quantity $X$ we define the flux surface average as

[TABLE]

and the transit average as

[TABLE]

with $d\tau=d\theta{B_{0}}/({v_{||}}\mathbf{B}_{0}\cdot\nabla\theta)$. We allow $X$ to depend on $v_{\|}$ and take the $\theta$ integrations in both averages (numerators and denominators) to be over a full poloidal circuit following a charged particle. In this way $v_{\|}$ and $\theta$ change signs together at a turning point for trapped particles, and odd functions of $v_{\|}$, such as $\overline{{v_{||}}}=0$ and $\overline{{v_{||}}/{B_{0}}}=0$, result in a vanishing transit average. Using $R/{R_{0}}(\psi)\simeq 1+\epsilon(\psi)\cos\theta$, $\epsilon=\epsilon(\psi)=r(\psi)/{R_{0}}$, $\underline{B}_{0}\equiv\left\langle{B_{0}}\right\rangle$, ${B_{0}}/\underline{B}_{0}\simeq{R_{0}}/R\simeq 1-\epsilon\cos\theta+\mathcal{O}({\epsilon^{2}})$, $\lambda=2\mu\underline{B}_{0}/{v^{2}}$, and

[TABLE]

we have for trapped particles that

[TABLE]

and

[TABLE]

For the passing particles we see from (119) and (120) that

[TABLE]

and

[TABLE]

when $W(\psi,\theta)$. For trapped particles $\overline{\xi W}=0$, while $W=1$ gives the passing result

[TABLE]

We employ a Shafranov shifted circular flux surface model (Shafranov, 1966; Helander & Sigmar, 2005). We retain the Shafranov shift ${\Delta_{S}}$ by taking

[TABLE]

with ${\Delta^{\prime}_{S}}=d{\Delta_{S}}/dr\sim\epsilon$. Then ${B_{0}}/\mathbf{B}_{0}\cdot\nabla\theta\simeq q{R_{0}(\psi)}\left[1-{\Delta^{\prime}_{S}}\cos\theta+\mathcal{O}({\epsilon^{2}})\right]$.

There are many integrals that need to be performed. The simple ones involve combinations of powers of $v_{\perp}$ and $v_{\|}$ multiplied by a Maxwellian without a transit average. More complicated ones involve transit and/or flux surface averages, for example, integrals of the form $\smallint{d^{3}}v{f_{0}}{v_{||}}\overline{{v_{||}}}$. The integrals involving transit average are most conveniently performed using $v$ and $\lambda$ variables so that

[TABLE]

with $\sum\nolimits_{\sigma}{\sigma}|{v_{||}}|\to 2|{v_{||}}|$ to account for both directions of ${v_{||}}$ for integrals over even functions of ${v_{||}}$.

We also see from (119) and (120) that

[TABLE]

and

[TABLE]

where we use (124). The preceding gives, for example,

[TABLE]

In the preceding evaluations and hereafter we will often make use of

[TABLE]

and

[TABLE]

for arbitrary gyrophase independent functions $X$ and $Y$.

Rosenbluth & Hinton (1998) analytically evaluated $\left\langle\smallint{d^{3}}v{f_{0}}{v_{||}}({v_{||}}-\overline{{v_{||}}})\right\rangle$ to find

[TABLE]

where the numerical constant $\gamma\approx 1.64$ comes from $3\int_{0}^{1-\epsilon}{d\lambda\overline{\xi}}=2\left[1-\gamma{\epsilon^{3/2}}+\mathcal{O}({\epsilon^{2}})\right]$. It is instructive to obtain their result by first using

[TABLE]

where

[TABLE]

Then using (126) for the passing, with $\bar{\xi}=0$ for the trapped, and $\cos\theta=1-2{\sin^{2}}(\theta/2)$ gives

[TABLE]

We then let $\alpha=\theta/2$ and introduce

[TABLE]

with

[TABLE]

to obtain the passing result in terms of a complete elliptic integral of the first kind

[TABLE]

Using $d\lambda=4\epsilon kdk/{\left[(1-\epsilon){k^{2}}+2\epsilon\right]^{2}}$ then yields

[TABLE]

Inserting $\oint d\theta/\xi$ from (140) gives

[TABLE]

Using

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

gives the result needed to recover (134)

[TABLE]

Using (131) and (147) allows us to generalize the Rosenbluth and Hinton result to find

[TABLE]

This $\theta$ dependence can also be checked by using $v_{||}^{2}={v^{2}}(1-\lambda{B_{0}}/\underline{B}_{0})$ and

[TABLE]

to see that

[TABLE]

which when integrated agrees with (147).

We can account for the different behaviour of the trapped (and barely passing) and passing particles by using $\overline{{v_{||}}}=0$ and estimating ${v_{||}}\sim{\epsilon^{1/2}}v$ and ${\int d^{3}}v\propto{\epsilon^{1/2}}$ for the trapped particles, while using $\overline{{v_{||}}}\simeq{v_{||}}\left[1+\mathcal{O}(\epsilon)\right]$ for the passing ones. In (148) these estimates give the order $\epsilon$ poloidal variation as coming from the passing particles while the order $\epsilon^{3/2}$ behaviour is due to the trapped.

However, we have to be more careful with a related integral since the Shafranov shift will enter. Using

[TABLE]

for the trapped, while noting that the passing particle result depends on the Shafranov shift we find

[TABLE]

where we have used (126). Consequently, when integrated over $\lambda$ we obtain

[TABLE]

Notice that (153) is $\mathcal{O}(\epsilon)$ since our shifted circle model requires ${\Delta^{\prime}_{S}}=d{\Delta_{S}}/dr\sim\epsilon$. Consequently, we see that $\overline{\xi\cos\theta}\sim\epsilon$ for the passing particles as expected.

A related integral involves $\underline{B}_{0}^{-1}\overline{{B_{0}}\xi}$, with $\underline{B}_{0}^{-1}\overline{{B_{0}}\xi}=0$ for the trapped, and

[TABLE]

for the passing particles, where we can use (126) and (140). As a result, expanding and neglecting ${\epsilon^{2}}$ corrections yields

[TABLE]

where we use $\overline{\cos\theta}\sim\epsilon$ for the passing and $\overline{\cos\theta}\sim 1$ for the trapped particles.

Another integral of interest is

[TABLE]

with $\bar{\xi}=0$ for the trapped particles, $\bar{\xi}$ given by (126) and (140) for the passing particles, and

[TABLE]

Hence, we need

[TABLE]

along with

[TABLE]

and

[TABLE]

Therefore, we find

[TABLE]

As a result,

[TABLE]

To perform some of the more complicated integrals we need to evaluate some transit averages. We have already evaluated $\oint{d\theta}{\xi^{-1}}$ for the passing particles, but now we also need it for the trapped particles. For the trapped we let

[TABLE]

with

[TABLE]

and

[TABLE]

then $\sqrt{{\kappa^{2}}-{{\sin}^{2}}(\theta/2)}=\kappa\cos\alpha$ and $\cos(\theta/2)(d\theta/2)=\kappa\cos\alpha d\alpha=\sqrt{1-{\kappa^{2}}{{\sin}^{2}}\alpha}d\theta/2$, give the half bounce result

[TABLE]

We also need to evaluate $\overline{\cos\theta}$ for the trapped particles, but we also evaluate it for the passing particles to check our estimates. For the passing we use ${k^{2}}\cos\theta={k^{2}}-2+{\kern 1.0pt}2\left[1-{k^{2}}{\sin^{2}}(\theta/2)\right]$ to find

[TABLE]

where $E$ is a complete elliptic integral of the second kind. For the freely passing particles $\lambda\to 0$ giving ${k^{2}}\to 2\epsilon\lambda\to 0$ so we recover the estimate $\overline{\cos\theta}\sim\epsilon$. For the trapped particles we again use (163)–(165) to find

[TABLE]

where we recover the estimate $\overline{\cos\theta}\sim 1$. The preceding gives the passing particle result

[TABLE]

and the half-bounce trapped particle result

[TABLE]

In addition, we will need $\oint{d\theta{\xi^{-1}}{{\sin}^{2}}\theta}$ for the trapped and passing. Using

[TABLE]

gives the half-bounce trapped particle result

[TABLE]

where we use #2.583.4 on p. 182 of Gradshteyn & Ryzhik (2007). We will also need to evaluate $\oint{d\theta{\xi^{-1}}{{\sin}^{2}}\theta}$ for the passing (because of the barely passing particle contribution). Using

[TABLE]

gives the passing particle result

[TABLE]

Using the preceding we can evaluate more complicated integrals like

[TABLE]

where we note that

[TABLE]

To evaluate any $\mathcal{O}({\epsilon^{3/2}})$ terms we Fourier decompose keeping only the leading harmonic,

[TABLE]

Then we can determine $\varsigma$ from

[TABLE]

Using the preceding results for $\oint d\theta{\xi^{-1}}\cos\theta$ and $\overline{\cos\theta}$ yields

[TABLE]

where the last integral is well behaved since ${\left[2E(k)-(2-{k^{2}})K(k)\right]^{2}}\propto{k^{8}}$ at small $k$, and the $k$ integral is the passing particle contribution and the $\kappa$ integral the trapped particle contribution. Consequently,

[TABLE]

Combining (162) and (178) gives

[TABLE]

Another related integral that we will need is

[TABLE]

Using $(v_{||}^{2}+\overline{v_{||}^{2}})-2{v_{||}}\overline{{v_{||}}}=2(\overline{v_{||}^{2}}-{v_{||}}\overline{{v_{||}}})+(v_{||}^{2}-\overline{v_{||}^{2}})$ and recalling that to evaluate (184) we used

[TABLE]

then we see that

[TABLE]

We also require the $\mathcal{O}(\epsilon)$ flux function

[TABLE]

where we have used (136), (152), and (161).

The more complicated, poloidally dependent integral $\smallint{d^{3}}v{f_{0}}v_{\bot}^{2}\overline{v_{||}^{2}\cos\theta}$ is required to $\mathcal{O}(\epsilon^{3/2})$. Flux surface averaging gives

[TABLE]

We also need the $\cos\theta$ dependent $\mathcal{O}(\epsilon^{3/2})$ portion of this poloidally dependent integral. Fourier decomposing by keeping only the fundamental

[TABLE]

we evaluate the Fourier coefficient by multiplying by $q{R_{0}}\underline{B}_{0}^{-1}\cos\theta\mathbf{B}_{0}\cdot\nabla\theta$ and flux surface averaging to find

[TABLE]

Inserting $\oint{d\theta{\xi^{-1}}{{\sin}^{2}}\theta}$ and $\overline{\cos\theta}$ for the trapped and passing particles gives

[TABLE]

where in the last integral $\left[2E(k)-(2-{k^{2}})K(k)\right]\left[(2-{k^{2}})E(k)-2(1-{k^{2}})K(k)\right]\propto{k^{8}}$ at small $k$ to keep it well behaved, and the $k$ integral is the passing particle contribution and the $\kappa$ integral is the trapped particle contribution (a notation used from here on). Then the preceding gives

[TABLE]

We next define the more involved integral I and approximate it by

[TABLE]

with $\sigma$ and $\upsilon$ order unity constants. The form $I$ is then rewritten using

[TABLE]

and

[TABLE]

to find that

[TABLE]

and

[TABLE]

with the domain of integration of order $\epsilon^{1/2}$ for the trapped particles and order unity for the passing particles. We use $\overline{{v_{||}}}=0$ for the trapped particles, estimate ${v_{||}}\sim{\epsilon^{1/2}}v$ for the trapped (and barely passing) particles, and use $\overline{{v_{||}}}\simeq{v_{||}}\left[1+\mathcal{O}(\epsilon)\right]$ for the freely passing particles. Therefore, we anticipate that the trapped particle ($\sim\epsilon^{5/2}$) contributions will give the final form of (194), with passing particle contributions $\sim\epsilon^{3}$.

To verify this more completely we next form

[TABLE]

where $\theta$ derivatives of the $\lambda$ limits do not contribute because terms containing transit averages have $\theta$ independent $\lambda$ limits and those not containing transit averages are multiples of ${v_{||}}$ that will vanish at the upper limit of $\lambda=\underline{B}_{0}/{B_{0}}$. Using $v_{||}^{2}={v^{2}}(1-\lambda{B_{0}}/\underline{B}_{0})$ and ${v_{||}}\partial{v_{||}}/{\kern 1.0pt}\partial\theta=-({v^{2}}\lambda/2\underline{B}_{0})\partial{B_{0}}/\partial\theta\simeq-\epsilon({v^{2}}\lambda/2)\sin\theta$ gives

[TABLE]

where we have used (162) and (184). Integrating gives to lowest order

[TABLE]

giving $\upsilon=5\gamma$ and $\sigma$ given by (199) to be the constant found by evaluating

[TABLE]

The $\varsigma\cos 2\theta$ term in (202) is of no consequent for our purposes since we only keep $\cos\theta$ terms.

We also need

[TABLE]

where we define $h=\underline{B}_{0}/{B_{0}}\simeq 1+\epsilon\cos\theta$ and note that for the trapped particles $\overline{{v_{||}}h}=0$. Using ${v_{||}}h-\overline{{v_{||}}h}\simeq{v_{||}}-\overline{{v_{||}}}+\epsilon{v_{||}}\cos\theta-\epsilon\overline{{v_{||}}\cos\theta}\simeq{v_{||}}-\overline{{v_{||}}}+\epsilon{v_{||}}\cos\theta+\mathcal{O}({\epsilon^{2}})$ we expand for $\epsilon\ll 1$ to find

[TABLE]

where we use ${v_{||}}-\overline{{v_{||}}}\sim v\epsilon$ for the passing and ${v_{||}}\sim v{\epsilon^{1/2}}$ for the trapped particles.

Finally, we can simplify (203) further by first forming $\oint{d\theta\xi}$. For the passing particles we let $\alpha=\theta/2$ to find

[TABLE]

while for the trapped we let $\sin(\theta/2)=\kappa\sin\alpha$ to obtain the half bounce result

[TABLE]

where we use

[TABLE]

from no. 2.584.6 on p. 186 of Gradshteyn & Ryzhik (2007) or use our previous results.

Inserting (140), (166), (206) and (207) into (203) and using $\bar{\xi}=0$ for the trapped and (126) for the passing particles gives

[TABLE]

Simplifying by noting that $\overline{{{(\xi-\bar{\xi})}^{2}}}\propto{k^{2}}$,

[TABLE]

We can simplify the trapped particle contributions further by using no. 2.584.15 on p. 187 of Gradshteyn & Ryzhik (2007),

[TABLE]

to write

[TABLE]

Then the first trapped particle term becomes

[TABLE]

upon using $\int_{0}^{1}d\kappa\kappa E(\kappa)=\frac{2}{3}$, $\int_{0}^{1}d\kappa{\kappa^{3}}E(\kappa)=\frac{{14}}{{45}}$, $\int_{0}^{1}d\kappa\kappa K(\kappa)=1$, $\int_{0}^{1}d\kappa{\kappa^{3}}K(\kappa)=\frac{5}{9}$, and $\int_{0}^{1}d\kappa{\kappa^{5}}K(\kappa)=\frac{{89}}{{225}}$, from pp. 615-616, no. 5.112.3-7 of Gradshteyn & Ryzhik (2007). To simplify the second trapped particle term we use

[TABLE]

to find

[TABLE]

where we use the preceding results to perform some of the integrals.

We can also simplify the simpler of the passing particle terms using

[TABLE]

and

[TABLE]

giving

[TABLE]

The remaining passing particle term can be rewritten as

[TABLE]

by using the passing particle results

[TABLE]

[TABLE]

and

[TABLE]

Putting all this together (210) becomes

[TABLE]

Numerically evaluating the integrals yields the constant $\sigma$ to be

[TABLE]