Multiscale Gyrokinetics for Rotating Tokamak Plasmas II: Reduced Models for Electron Dynamics

In this paper, we extend the multiscale approach developed in [Abel et. al., Rep. Prog. Phys., submitted] by exploiting the scale separation between ions and the electrons. The gyrokinetic equation is expanded in powers of the electron to ion mass ratio, which provides a rigorous method for deriving the reduced electron model. We prove that ion-scale electromagnetic turbulence cannot change the magnetic topology, and argue that to lowest order the magnetic field lies on fluctuating flux surfaces. These flux surfaces are used to construct magnetic coordinates, and in these coordinates a closed system of equations for the electron response to ion-scale turbulence is derived. All fast electron timescales have been eliminated from these equations. We also use these magnetic surfaces to construct transport equations for electrons and for electron heat in terms of the reduced electron model.


Introduction
The energy, particle, and momentum confinement of present-day fusion experiments and proposed future devices is limited by turbulent, rather than collisional transport. The turbulence that causes this transport has an essentially multiscale character. It occurs on (perpendicular) spatial scales smaller than those associated with the equilibrium and on timescales shorter than those associated with transport of mean quantities but much longer than the gyroperiod. This is the basis for the multiscale gyrokinetic approach presented in [1,2], characterised by the small parameter = ρ/a, where ρ is the thermal gyroradius and a is a typical equilibrium length scale. In these treatments, all species are considered equal. However, there is a second scale separation that is important, the scale separation between the plasma ions and the electrons, characterised by the small parameter δ = m e /m i , where m e and m i are the electron and ion masses respectively. For a deuterium plasma δ ≈ 1/60.
Ion-scale fluctuations (e.g. Ion Temperature Gradient (ITG) [3,4,5], Parallel Velocity Gradient (PVG) [6,7] or Trapped Electron Mode (TEM) [8] driven turbulence) have length scales perpendicular to the magnetic field comparable to the ion gyroradius and have frequencies comparable to the ion transit or bounce frequencies. In contrast, electron-scale fluctuations have perpendicular length scales comparable to the electron gyroradius, and frequencies comparable to the electron transit or bounce frequency. Unless they are suppressed (e.g. by strongly sheared flows), ion-scale fluctuations with their larger "mixing length" usually dominate the transport -and are therefore considered to be more important. In this paper we focus on the electron response to such fluctuations.
The electrons are much lighter than the ions, hence ion-scale fluctuations are longwavelength compared to the electron gyroradius and low-frequency when compared to the electron transit or bounce frequencies. It is clear that in simulating ionscale turbulence one does not wish to resolve the short electron timescales. This has motivated the development of models for the electron response to ion-scale turbulence that somehow eliminate fast processes due to electrons [9,10,11,12]. In this paper, we derive for the first time a set of equations governing the electron behaviour in the presence of ion-scale fluctuations. These equations eliminate the fast electron timescale via a rigorous expansion of the gyrokinetic equations in δ 1, subsidiary to the expansion in 1. The fast processes arise from the streaming of electrons along magnetic field lines. In electrostatic turbulence simulations, this streaming is along the unperturbed mean field and easily handled. The electrostatic approximation only holds at very low β (the ratio of thermal to magnetic pressure, and a measure of the importance of magnetic perturbations). However, any future reactor design will operate with as high a β as possible. Indeed, spherical tokamaks (such as MAST [13]) already operate at finite β. Thus, we require an electron response model that incorporates the fluctuations of the magnetic field in a fully general way. As we shall see, if we consider only ion-scale 1.5-1.6% 1.5-1.6% 1.5-1.6% Table 1. Typical length and time scales in selected fusion devices.
fluctuations then the topology of the magnetic field is preserved by the fluctuations (see Section 3.2). However, to determine the electron behaviour we will need to make assumptions about the initial structure of the magnetic field and the effect of other scales. In particular, whether the field is regular, with nested flux surfaces, or stochastic. We show in Section 3.3, that if the magnetic field were significantly stochastic (i.e. if the ion heat transport due to the stochastic field were comparable to the transport due to the fluctuating E × B drifts) then the electron heat transport would be larger than the ion transport by a factor of δ −1 1. Such poor electron heat confinement is not observed in experiments, indeed it is common for the level of ion heat transport to exceed the level of electron heat transport in conventional tokamaks. Thus, we assume that any stochastic component of the magnetic field is small enough to be ignored (see Section 3.3) and therefore that the field lines lie on fluctuating magnetic surfaces.
The structure of the rest of the paper is as follows. In Section 2.1, we introduce the electron gyrokinetic equation and attendant notation. In Section 2.2, we develop the formal expansion in δ 1 and order all physical quantities with respect to this small parameter. The key result of this paper -the gyrokinetic system with reduced electron equations -is derived via a systematic expansion in δ. This derivation is performed in Sections 3 and 4 with the results concisely summarised in Section 5.
In Section 3.2, we use the expansion in δ 1 to prove that long-wavelength lowfrequency fluctuations preserve the topology of magnetic field lines, generalising the results of [14,15,16,17]. With the assumption of an initially regular magnetic field, this result enables us to introduce a coordinate system aligned to the exact magnetic field in Section 3. It is this set of coordinates that allows us, in Sections 4.2 and 4.3, to average over the fast electron motion along exact magnetic field lines and so close our equations for the electron distribution function. These equations are particularly simple in the limit of strong collisions (ν e ω) which we derive in Section 6. The second major result of this paper is a simplified system of transport equations for electrons. This is presented in Section 7 and derived in Appendix H. These equations emphasise the particularly simple form that the transport fluxes take in the low-massratio limit. We also present the transport equations for the collisional electron model of Section 6, which take on an even simpler form.
We draw the results of the paper together in Section 8. We make a detailed comparison of our equations to previously-derived models in Section 8.1, emphasising both the differences and the similarities. Finally, we conclude in Section 8.2.

Electron Gyrokinetics and the Mass-Ratio Expansion
We study the response of electrons to ion-scale fluctuations within the framework of multiscale gyrokinetics [1]. First, we will introduce this framework in Section 2.1 along with all attendant notation. Then, in Section 2.2, we formally order all quantities in the small parameter δ = m e /m i .

Gyrokinetics in a Rotating Tokamak Plasma
In this section, we introduce the assumptions and resulting equations of this formalism. The detailed derivation of these results can be found in [1].
First, we split all physical quantities into mean and fluctuating parts: where B is the magnetic field, E the electric field, f s the distribution function of species s and · turb is some average over the fluctuations. We assume that the fluctuations obey the standard gyrokinetic ordering ‡ : where ω is the typical fluctuation frequency in the frame rotating with the plasma, Ω s = Z s eB/m s c is the cyclotron frequency of species s with charge Z s e and mass m s in the mean magnetic field B = |B|, k and k ⊥ are typical wavenumbers parallel and perpendicular to the mean magnetic field, ν s is a typical collision frequency, ρ s = v ths /Ω s is the thermal Larmor radius, v ths = 2T s /m s is the thermal speed of species s with mean temperature T s and a is a typical equilibrium length scale. The mean quantities (B, E, F s ) are assumed to vary on the long length scale a and the long (transport) timescale τ E , the energy confinement time. This long timescale is ordered so that ) ‡ In this paper, as in [1], the symbol ∼ is used to mean "is the same order as" (with respect to the appropriate expansion) rather than the more usual "is asymptotically equivalent to".
Typical values for the timescales and length scales in a tokamak are given in Table (1), which shows the very large scale separation, characterised by the small parameter . This scale separation in time and space of the mean and fluctuating quantities motivates us to define the average · turb in terms of temporal and spatial averages over intermediate time and length scales. This is done precisely in equations (14)- (18) of [1].
It is shown in Section 3.3 of [1] that if we introduce the cylindrical coordinate system (R, z, φ) and assume axisymmetry with respect to φ, the mean magnetic field can be written (to lowest order in ) as where the poloidal flux function ψ is given in terms of the mean vector potential A by and I(ψ, t) = R 2 B · ∇φ is the toroidal component of the mean magnetic field. We will also need the following alternate form of the magnetic field [18,19]: where α is a variable that labels different field lines within a flux surface. We cut the toroidal surface on the inboard side to force α to be single-valued. This cut forms an axisymmetric ribbon with one edge the magnetic axis and the other the plasma boundary or separatrix. Letting l be the distance along a field line, (ψ, α, l) is a set of field-aligned coordinates for B. There are many different forms of field aligned coordinates and none of our results depend on a particular choice. Nevertheless, we will use (ψ, α, l) for specificity. Axisymmetric functions are functions of ψ and l but not α. Note that both ψ(R, z) and I(ψ) are mean quantities and therefore vary on the transport timescale. We assume that the mean electric field is large, E ∼ (v ths /c) B, which can be shown to imply that the plasma rotates toroidally with a velocity that is species-independent and whose angular velocity only depends on the flux label ψ [20] and the slow transport timescale: If we express the mean electric field in terms of potentials using Gaussian units and Coulomb gauge, then the scalar potential takes the form and we have the following relationship between ω(ψ, t) and Φ(ψ, t): With these results for the fields, our orderings imply that the distribution function f s of species s takes the following form: Let us explain this notation. We have followed [1] and changed phase-space variables from (r, v) to (R s , ε s , µ s , ϑ, σ): where b = B/B is the unit vector along the mean magnetic field, the peculiar velocity w is where e 1 and e 2 are two arbitrary orthogonal unit vectors perpendicular to the magnetic field, the direction of the parallel motion is σ = w /|w |, and we have defined The quantity ψ * s (r, v, t) is equal to c/(Z s e) times the canonical toroidal angular momentum and is thus conserved in an exactly axisymmetric system -i.e., in the absence of fluctuations. In these variables, (14) splits the mean distribution function into the near Maxwellian equilibrium and the neoclassical distribution function F 1s , which describes large-scale O( ) deviations from a Maxwellian. The function N s is related to the mean density n s by We now drop the explicit slow time dependence of all quantities -this is for notational clarity only and is not an assumption that the mean quantities are constant. The fluctuating distribution function, given by (15), is composed of the Boltzmann response (the first term of (15)), where the fluctuating potentials δϕ and δA have been combined to form the electrostatic potential in the frame rotating with the plasma and the gyrokinetic distribution function h s . Assuming that all mean quantities are known, then the gyrokinetic distribution function h s is given by the gyrokinetic equation [1,2]: where the guiding-centre drift velocity is and where we have defined the gyrokinetic potential and the fluctuating velocity due to χ § We have also introduced the gyroaverage at constant R s , ε s and µ s : The equation for h s is closed through Maxwell's equations for the fluctuating fields. The fluctuating potential δϕ obeys the quasineutrality condition: where the integral over velocities is performed at constant r and the gyroaverage at constant r, w , and w ⊥ is defined by The fluctuating magnetic field is determined from δA = b · δA and δB = b · δB, which obey the parallel component of Ampère's law and the perpendicular part of Ampère's law (see the discussion in Section 7.5 of [1]) respectively. The evolution equation (24) and constituent relations (29), (31) and (32) form a closed system of equations for determining h s , δϕ , δA , δB on the turbulent timescale. Typically, the system (24), (29), (31), and (32) is solved until statistical equilibrium is reached (a few nonlinear turnover times). The slowly-varying equilibrium is treated as constant during this evolution. The slow time dependence of F 0s is determined by transport equations for ω(ψ, t), N s (ψ, t), T s (ψ, t) and I(ψ, t) and the equilibrium condition determining ψ(r, t). Such equations are given in [1].

The Mass-Ratio Expansion
In order to express the scale separation of the electron and ion scales, we introduce the secondary small parameter δ = m e /m i 1. We will use a subscript i to denote any ion species, of which there may be many. We consider the expansion in δ 1 to be a subsidiary expansion of the gyrokinetic system (equations (24), (29), (31), and (32)) . Thus, we assume that = ρ i /a, i.e., that the fundamental gyrokinetic expansion parameter is that of the ions, and that All other dimensionless parameters are treated as finite -i.e., independent of δ. Therefore, we are assuming that T i ∼ T e and β = 8πp/B 2 ∼ 1, where p is the plasma pressure. In typical fusion plasmas, β is often small and one might worry about this assumption. A more detailed analysis shows that our expansion in δ requires that [21] β m e m i = δ 2 , which is equivalent to requiring that the Alfvén speed is smaller than the electron thermal speed. ¶ In plasmas of interest, this is indeed satisfied and we are justified in treating β as finite. We assume that fluctuating quantities vary on ion rather than electron scales. Thus, we order the timescales of fluctuating quantities in the rotating frame and the size of the fluctuations by ∂ ∂t so that the typical fluctuation timescale, the ion parallel streaming time, and the ion nonlinear timescale due to the fluctuating fields are all comparable. We assume that the toroidal rotation velocity u is comparable to the ion thermal speed u ∼ v th i , which implies that the electric field is ordered as This is in contrast to the approach taken by [12], see Section 8.1 for details. ¶ For cases where β is so low that this does not hold, an expansion of slab gyrokinetics for very low β plasmas has been carried out in [22].
The perpendicular length scales are of course ordered as Clearly, the turbulent dynamics of the ions are not simplified by the subsidiary expansion in δ and hence h i is determined by the gyrokinetic equation (24). The electron dynamics, however, are simplified because so the parallel electron motion is much faster than the evolution of the fluctuations + and also so we can neglect finite-electron-gyroradius effects. We choose the electron collision rate ν e to be comparable to the fluctuation frequency, i.e., Note that although electron collisions are comparable to the fluctuation frequency, and, therefore, important in the electron fluctuation dynamics, this would conventionally be considered a low collisionality plasma since ν eff = ν i a/v th i ∼ ν e a/v the ∼ δ 1. Under these assumptions, we expand h e , δϕ , δA and δB as h e = h (0) e + h (1) e + · · · , δϕ = δϕ (0) + δϕ (1) + · · · , δA = δA (0) + δA (1) + · · · , and δB = δB (0) + δB (1) where h

Zeroth Order: Particular Solutions, Flux Conservation and Field-Aligned Coordinates
In this section, we apply the orderings of Section 2.2 to the gyrokinetic equation for electrons and examine the consequences to lowest order.
In order to expand the gyrokinetic equation, in Appendix A, we express χ R explicitly in terms of the fields and find, (A.3): + Examining the assumption that ω/k v the ∼ δ, we might worry that if k is small then this assumption would be violated. Thankfully, studies of strongly nonlinear turbulence (both numerically-simulated ITG and experimentally measured turbulence) suggest that it is, in fact, critically balanced [23,24] all timescales are comparable and the nonlinear processes fix k such that ω ∼ k v thi . Thus, we do not expect our assumption to break down in practice. * In fact, we will only require δϕ or δB to lowest order in δ. Thus, we will drop the superscripts on these fields.
Substituting this result and the expansion (41) into the gyrokinetic equation (24), we find, to lowest order in δ: where δA is evaluated at R e and we have used u(R e ) · ∂w /∂R e = 0 (axisymmetry). These equations are written in terms of the variables R e , ε e , µ e and ϑ. It will be convenient to instead convert back to using r as our spatial variable. Since (43) is only accurate up to corrections of order O(δ 2 Ω i F 0e ) and all functions in (43) are evaluated at R e , we can use the fact that k ⊥ ρ e ∼ δ to replace R e by r everywhere. The function ε e is now given by the simplified expression: Defining b = B/| B| to be the unit vector along the exact magnetic field, we have Thus, after a little algebra (using (45)), we can rewrite all but one term of (43) as derivatives along the exact field line: where we have used the Maxwellian form (21) of F 0e , and divided through by F 0e . Equation (46) is an inhomogeneous ordinary differential equation along the exact field lines for δf e . In the remainder of this section, we will find the general solution of (46). The particular solution to the inhomogeneous problem is constructed in Section 3.1. In order to solve the homogenous part of (46), it is necessary to discover the structure of the exact magnetic field. Thus, in Section 3.2, we prove that magnetic flux is conserved to lowest order in δ and hence that the topology of the magnetic field is fixed. This result then allows us in Section 3.4 to introduce field-aligned coordinates in which it is simple to solve the homogenous part of (46). This general solution will provide us with an equation for δA and place constraints upon the form of h

Particular Solutions of (46)
We can now find particular solutions of (46). To do this, we write δf e in the following form: δf e = λ(r, ε e , µ e , σ)F 0e + g e (r, ε e , µ e , σ) where g e satisfies the homogenous equation This can clearly be done for any δf e by suitable choice of λ. Substituting (48) into (46) we obtain To find the form of λ, we take the derivative of this equation with respect to µ e and the second derivative with respect to ε e to obtain b · ∇ ∂λ ∂µ e = 0 and b · ∇ ∂ 2 λ ∂ε e 2 = 0, respectively. Clearly, any solution to these equations is either independent of µ e and a linear function of ε e or satisfies b · ∇λ = 0. Any solution of the second kind can be absorbed into the function g e , leaving us with the general solution for λ given by where we have written the linear function of ε e in the natural way as a perturbed Maxwellian. However, it should be noted that δn e (r) is not the perturbed electron density and δT e (r) is not the perturbed electron temperature because g e can have both density and energy moments. Inserting (52) into (50) and taking a derivative with respect to ε e , we find the following equation for δT e : b · ∇ ln T e + δT e T e = 0.
Finally, substituting both (52) and (53) into (50), we obtain: Familiar physics is contained in these two equations. It will be shown in Section 3.5 that the contribution to g e from passing electrons is a function only ofψ and not of α orl, and, therefore, so are its density and temperature. Thus, for the purposes of interpreting (54) and (53), we can add the density and temperature of the passing part of g e to δn e and δT e respectively. Hence, (53) describes how (part of) the temperature of the electrons is equalised along field lines, due to rapid thermal conduction. The only part of the temperature that does not equalise along the field line, as we shall see when we solve (49), is the part of the temperature of the trapped particles contained in g e . This is to be expected; the trapped particles do not traverse the entire field line and so cannot equalise their temperature along it. Equation (54) describes how the pressure gradient of the electrons (again, not including the trapped particle pressure) along the field line is balanced by a parallel electric field. As we shall demonstrate, it is useful to think of (54) as an evolution equation for δA where δϕ and δn e are given.

Flux Conservation
To solve (49), (53) and (54) we will need to know the spatial structure of the total magnetic field. As (54) determines δA , it determines δB ⊥ and thence b. In this section, we will prove that (54) implies that magnetic flux is conserved and consequently that the evolution of the magnetic perturbation cannot change the structure of the magnetic field lines.
Magnetic flux is said to be conserved if there exists an effective velocity field u eff such that closed curves moving with this velocity always enclose the same amount of magnetic flux. It is proved by Newcomb in [25] that if such a u eff exists, it also preserves magnetic field lines and thus the magnetic topology is fixed. Now, if the parallel electric field satisfies where ξ is a single-valued scalar function, then the velocity field (where B = | B| and U (r) is an arbitrary single-valued function) satisfies E + u eff × B/c = −∇ξ. Thus Faraday's law becomes ∂ B/∂t = ∇ × (u eff × B). The familiar proof of the frozen-in theorem (that is reproduced in almost all MHD textbooks) clearly applies if we recognise u eff as the flux-and field-line-preserving velocity. The condition (55) is a special case of the general conditions for the existence of frozen flux that were investigated in [25].
In Appendix B, we prove that we can write E · b in terms of potentials as In Appendix C, we prove that the mean magnetic flux is conserved irrespective of the fluctuations, and so (see (C.7)) where n 1e is the O( ) correction to the mean electron density (see discussion before (C.7)). Finally, using (54) and (58) in (57), we obtain which is indeed in the form (55). † Thus, we can define ξ = T e e ln N e + δn e + n 1e n e + ϕ 0 +ξ, whereξ is an arbitrary single-valued function that satisfies and with this definition, (56) defines the field-line-preserving velocity u eff . In (56), we are free to choose U (r) andξ subject to the constraint b · ∇ξ = 0. We postpone making choices for these parameters until Section 3.5. Since the field is frozen to u eff the ionscale turbulence cannot alter the magnetic field structure on the turbulent timescale.

Structure of the Magnetic Field
In the previous section we have demonstrated that ion-scale turbulence preserves field lines. The topology of the field is thus fixed. But what is this topology? Electronscale fluctuations can alter the topology and so can generate a stochastic field ‡ from an initially regular one. Such fluctuations may arise "locally" from electron-scale instabilities (e.g. Electron Temperature Gradient instabilities or Microtearing), but electron-scale fluctuations may also be driven "non-locally" by a cascade from ion-scale fluctuations. In both cases, some tangling of the magnetic field is possible, indeed probable. Thus, at some intuitive level, the magnetic field is expected to be stochastic. A stochastic field will be advected by u eff and remain stochastic. If the magnetic field were stochastic and stationary then the heat diffusivity of the electrons would be given by the well-known Rechester-Rosenbluth formula [26]: where D st is a species-independent diffusion coefficient characterising the stochasticity of the field § and l c is the parallel correlation length of the stochastic component of the † Strictly speaking, we have proven that Substituting this into Faraday's law, we find that Estimating the size of the left hand side and the second term on the right hand side, we see that the correctionẼ will affect neither the lowest-order mean field B nor the lowest-order fluctuation δB. Thus, we can safely ignore it. ‡ A stochastic field is one in which points on neighbouring field lines become exponentially separated as they travel along field lines [26]. § The field-line diffusivity D st is defined as [26] where δr(l) is the radial displacement of a point after following a field line for a distance l and the overbar denotes statistical averaging over field lines (the factor of 2 in the denominator is chosen to agree with the definition in [26]).
field δB stoch . Let us assume the stochastic component of the field is comparable to the ion-scale field perturbations (i.e., δB stoch /B ∼ ) and the correlation length is the system size (i.e., l c ∼ a). Then, (62) predicts that the electron heat diffusivity in a tokamak should be δ −1 larger than the ion heat diffusivity (both the electrostatic and electromagnetic parts of the ion heat diffusivity). This is not backed up by experimental observations; ion and electron heat diffusivities are observed to be of the same order of magnitude. Thus, we take this as experimental justification for assuming that the amount of magnetic field stochasticity is limited. We therefore assume good flux surfaces in solving Eqs. (49), (53) and (54) -these surfaces will be approximations to the exact field that we will refer to as the "ion-scale field". The ion-scale turbulence will distort the ion-scale field but the surfaces will remain intact. A smaller stochastic field component with l c ∼ a and δB stoch /B ∼ δ 1/2 would cause observable levels of transport and we cannot rule out a component of this size. Indeed, such a field would change the final electron equations in this paper. However, we will assume for simplicity that the stochastic component of the field is zero to order O( δ) (i.e., δB stoch /B δ) in this paper.

Field-Aligned Coordinates
If we write the ion-scale magnetic field in the Clebsch form thenψ is a single-valued function that labels the flux surfaces andα is a function that labels different field lines within a given flux surface. To ensure thatα is single-valued, we make a cut in the poloidal domain on the inboard side of the tokamak. This cut makes a ribbon forming a toroidal loop with one edge of the ribbon being the magnetic axis and the other the plasma boundary or separatrix. We allow the cut to move with the velocity u eff . This is legitimate as u eff is continuous across the cut. Since the cut is perturbed by the turbulence, it will not now be axisymmetric. The expression (63) for the magnetic field allows us (see, e.g., [19]) to introduce a system of field-aligned coordinates (ψ,α,l) wherel measures distance along a given field line. By definition, l runs from one side of the cut in the poloidal plane to the other side of the cut. As We can formulate this condition precisely as a condition on the field-line diffusivity (defined in the footnote on page 13) . The limit on stochasticity can be expressed as where we obtain this estimate by insisting that the change in a function g following the field line is dominated by explicit variation along a field line (k g) rather than the field line connecting neighbouring points in the perpendicular plane (D st k 2 ⊥ g). The existence of a stochastic field giving rise to a field-line diffusivity of size D st ρ e would not invalidate the existence ofψ andα, as (for the purposes of studying ion-scale turbulence) we could defineψ surfaces by the mean radial location of a field line averaged along a magnetic correlation length l c . This introduces a displacement ρ e of theψ surface from the exact position of the field line, but for turbulence where k ⊥ ρ e 1 this correction is negligible. δB B, the variables aligned to the exact field only deviate slightly from those aligned to the background field: Importantly, this means that all axisymmetric mean quantities are, to lowest order in , functions ofψ andl but not ofα. We need expressions for the spatial derivative operators b·∇ and the Poisson bracket {a, b} = b·(∇a × ∇b) in our new coordinates. By using the fact that b·∇ψ = b·∇α = 0, we discover that for any function a(r). Similarly, by using the fact that b · (∇α × ∇ψ) = B to lowest order, we see that the Poisson bracket takes the form Since u eff is a field-line-preserving velocity, we have that and Asl ∼ l and ∇ · u eff ∼ ω, u eff also approximately preserves the length of field lines and we have, to lowest order in , It is of course convenient to choose the cut in the definition ofα such that it is convected with u eff . Finally, the time derivative at constantψ,α andl is related to the time derivative at constant r by

Solution of the Lowest-Order Equations
We now use the magnetic coordinates introduced in Section 3.4 to solve (49) and (53). We split velocity space (ε e and µ e ) into two regions: the passing region where, for a givenψ, the equation has no solutions for anyl (i.e., ε e > µ e B − eϕ 0 for alll) and the trapped region in which there are two values ofl at which w = 0. The points where w vanishes are called the bounce points of the electron's orbit. In the passing region, (49) implies that ∂g e ∂l = 0, so g e is constant along field lines. As most magnetic surfaces are irrational, one field line spans the entire surface and so, for passing electrons, g e cannot depend onα: In the trapped region of velocity space, (49) implies that either w (l,ψ, ε e , µ e ) = 0, or ∂g e ∂l = 0.
Thus, g e must be constant along magnetic field lines, but only between the bounce points. So g e can depend on the location of the bounce points, and, therefore, onα: where the dependence on σ has been dropped as g e must be continuous as w passes through zero and changes sign. For convenience, we define g pe to be zero in the trapped region and g te to be zero in the passing region, so the complete solution is just the sum of the two functions: The solution of (53) in our new coordinates is for some functionT e . As we can absorb any part of δT e that is only a function ofψ into g e , we can pickT e (ψ) = T e (ψ) and so We can also now pick the arbitrary functionξ in the definition (60) of ξ to be any function ofψ. In Appendix E, we find convenient forms for the arbitrary functionsξ(ψ) and U (r) to obtain (see (E.6)) where we have defined With this definition, we can express the general solution of (46) as and the evolution equation (54) for δA becomes This equation allows us to solve for δA correctly to lowest order in δ. Equations for the evolution of ζ, g pe and g te at the next order in our expansion are obtained in Section 4.

Solution for δB ⊥
At this juncture, we have enough information to determine δB ⊥ in two ways. Firstly, we have the evolution equations (67) and (68) forψ andα, with u eff given by (79). This enables us to calculateψ andα correctly to first order in . In Appendix D, we demonstrate that this gives us enough information to construct δA from the equations Secondly, we could solve for δA from the evolution equation (82) and then use (83) to determineψ andα. Either of these methods is valid, and provides a solution for δA correct to the requisite order. For the purposes of simulating our equations, a third option presents itself: determine δA from (82),ψ andα from (67) and (68), and then use (83) as consistency checks upon the solutions thus obtained. We choose the first of these options, i.e., to eliminate δA in favour of the fluctuating fieldsψ − ψ andα − α.
One point remains to be resolved. It would appear that knowingψ andα determines δB completely (via (63)). However, to determine δB fromψ andα, one must knowψ andα to O( 2 ) or, equivalently, the compressible part of u eff to higher order (O( 2 v th i ); see Appendix D), but we do not. Therefore, we determine δB from (32) which involves only lowest-order quantities -quantities that we do calculate.

Parallel Electron Currents
As explained in the previous section, we have obtained a solution for δA . However, δA is usually determined from the parallel component of Ampère's Law, (31). Substituting our expansion for h e (41) into the parallel component of Ampère's law (31), we find where a sum over s = i denotes a sum over all ion species. In this equation, the term due to h (0) e is larger than the rest by δ −1 . Thus, to lowest order, This is a constraint on g pe in the lowest-order solution given by (81). Note that the trapped particles do not contribute to the current constraint. The parallel electron flow (and current) is comparable to that of the ions and contained in h Thus, instead of determining δA from the electron current, Ampère's law determines the electron flow from δA [21]. As we are no longer using δA to describe δB ⊥ in our system of equations, we consider ∇ 2 ⊥ δA to be a shorthand for a complicated expression involving derivatives ofψ − ψ andα − α via the relations (D.10) or (D.11). We require δu e to complete the solution in the next order, but we do not need any information about h (1) e other than δu e .

First Order: Dynamical Equations for Electrons
In this section, we return to the gyrokinetic equation (24) and proceed to the next order in the mass-ratio expansion in order to find evolution equations for g e and ζ.
Subtracting (43) from (24) and using the solutions for h (0) e and δA from Section 3.5, we obtain the first-order part of the gyrokinetic equation: where we have used (45), rewritten terms involving V χ R as Poisson brackets using (27) and (66), and, as in Section 3, replaced R e with r as our spatial variable. Thus, the collision operator in (88) is evaluated at R e = r. The drift velocity in (88) now contains only the leading-order part: We now rewrite (46) in terms of δf e and ζ: where we have used (47), (80) and the fact that which follows directly from (79) and (66).
Finally, we substitute the form of δf e from (81), to obtain where we have gathered all terms involving g e together on the left-hand side and all terms involving δA (1) on the right. Notationally, (92) can be confusing as we have made some abbreviations for the sake of conciseness. Firstly, all spatial derivatives are taken at constant t, ε e and µ e . Thus, for equations that are solved in the moving coordinate system As this notation is cumbersome, we will continue to use ∇ where appropriate (i.e., where it would expand to more than one derivative). Secondly, all the Poisson brackets in (92) are shorthand for the expression (66) in terms ofψ andα derivatives. Finally, this equation is, intrinsically, an equation in the frame moving with the magnetic field , and all functions should be interpreted as functions of t,ψ,α, andl rather than t and r. Equation (92) describes how g e , and thus δf e , evolves dynamically in response to given perturbed fields. As u eff is also considered known, (67)-(69) can be solved to give the complete transformation from the fixed to the moving coordinate system. Thus, we consider all the quantities in (92) to be known in both the fixed coordinates (which will be needed when we come to the field equations) (ψ, α, l), and the moving coordinates, (ψ,α,l). However, in its current form, (92) is not closed, as we know neither h (1) e nor δA (1) . In the following sections we derive evolution equations for ζ, g te and g pe as solubility constraints of (92).

Fluctuating Continuity Equation
The one piece of information that is known about h (1) e is its parallel velocity moment δu e , determined by (87). Thus, we integrate (92) over all velocities and find e carries no current, and We have also defined and We have also used ∇ · B = 0 and (D.11) to eliminate b and δA from our equation.
It is important to note that this equation is written in the moving coordinate system. It relates ζ(ψ,α,l, t) to g e and the other fluctuating fields, also given as functions of ψ,α andl. In this context,ψ − ψ should be viewed as just another fluctuating field. Therefore, all gradients should be thought of in terms of derivatives with respect toψ andα, for the Poisson bracket this is done via (66). That (94) is essentially the fluctuating continuity equation is not immediatly obvious. Let us attempt to make this more transparent. Clearly, the first two terms under the time derivative in (94) are the perturbed electron density due to the fluctuations; the time derivative of δB adds to this the fluctuation in the electron density due to the compression of the flux surfaces. ¶ This total density is changed ¶ In the derivation of the transport equations we show that so ∂δB /∂t is equal to the rate of perpendicular compression of the flux surfaces.
by various flows -the compression in the parallel flow δu e , the flow of electrons past the flux surface given by (c/B) b × ∇ζ, the grad-B drift due to δB , and the magnetic drifts. Indeed, even the final term in (94) can be interpreted as a compression -it is the compression of the electron fluid due to the perturbed flux surfaces "sticking out" into regions with different toroidal angular velocities. Equation (94) provides us with an evolution equation for ζ if we know g e and the perturbed fields. The evolution equations for g e given ζ are found in the next two sections.

Trapped Electrons and the Bounce Average
In order to derive an equation for g e in the trapped region of velocity space, we need to define an averaging operator that eliminates the fast electron bounce motion (i.e., annihilates w b · ∇ and thereby eliminates h (1) e ). The appropriate average is the bounce average along the perturbed field i.e., at fixedψ andα. For an arbitrary function a(ψ,α,l, ε e , µ e , σ) we define the bounce average to be where the integration is carried out from one bounce point where w = 0, along the magnetic field, to the second bounce point where w passes through zero and changes sign, and then back to the first bounce point. This integral is carried out at constant ψ,α, ε e , µ e and ϑ. Thus, equivalently, the bounce average is where l 1 and l 2 are the bounce points of the particle orbit. Clearly, the bounce average commutes with the time derivative at fixed moving coordinatesψ,α andl. In Appendix F, we prove that and that, to lowest order in , Using these results, we take (92) in the trapped region of velocity space and perform the bounce average to find where the last line in (92) has vanished under the bounce average because it is antisymmetric in σ (the sign of the parallel velocity). This equation is a closed equation for g te , once the fields and ζ are given. Note that all quantities must be expressed in the moving coordinates, in order to perform the bounce average. Again, as with (94), all gradients should be interpreted in terms of derivatives with respect to the moving coordinates as should the Poisson bracket (via (66)).

Passing Electrons and the Flux Surface Average
In the passing region of velocity space, the fast streaming of electrons along perturbed field lines causes a passing electron to sample an entire perturbed flux surface (ψ = const. surface). Thus, the averaging procedure we will need is the average over a perturbed flux surface. For any function a, this is defined by a(r, ε e , µ e , ϑ, σ) ψ = lim where the region of integration ∆ is the volume between the surface labelled byψ and that labelled byψ + ∆ψ, and the integral is taken at constant ε e , µ e and ϑ. In Appendix G, we prove that that the flux surface average commutes with the time derivative following the magnetic field, that and that, for any function a(r), We also prove that, for any fluctuating function λ, Taking (92) in the passing region of velocity space, multiplying by B/|w |, and flux-surface averaging (to eliminate w b · ∇h (1) e ), we obtain: ∂t where the final term in (92) vanishes because in the passing region g e is only a function ofψ and we can use (107). We still have to eliminate the final line of (109). We note that we only require the part of g pe that is even in σ -we require its density in the quasineutrality condition (29) and its pressure in perpendicular Ampère's law (32). Thus, we average (109) over the two signs of σ and henceforth g pe denotes only the even part. Hence we obtain This completes our solution for g e , as (110) is a closed equation for the part of g pe that is even in σ. It is useful to note that only the zonal (i.e.,α-independent) parts of ζ and δB will give a significant drive for g pe -theα-dependent parts will tend to cancel under the flux-surface average.

Field Equations
In this section, we finally close our system of equations by writing the field equations in terms of our solution for δf e and ion quantities. We find δϕ by using the solution (81) for δf e in the quasineutrality condition (29): Similarly, we find δB by using (81) in (32) to find These two equations determine δϕ and δB .
It is instructive to note that (111) and (112) are not naturally solved in the moving coordinate system. It is because of this coupling to the fixed coordinates (and the ion physics) that we require the mapping between the coordinate systems (from the solution of (67)-(69)).

Summary of the Reduced Electron Model
The set of electron equations we have derived is complicated: the electron time-scale has been removed at the expense of introducing more dependent variables and more equations, some integro-differential. In these equations, we have the nine unknownsψ, α, δu e , ζ, g te , g pe , h i , δϕ , and δB . We have provided a complete set of nine equations to determine these quantities. These are: • (67) and (68) to evolveψ andα which define the moving coordinate system that moves with velocity u eff given by (79), • (87) to determine δu e , • and the coupled set (94), (103), (110), and (24) to evolve ζ, g te , g pe , and h i with the constituent relations (111) and (112) determining δϕ and δB .
Together, these equations comprise a rigorously-derived model for ion-scale turbulence, accurate up to corrections of order δ. The only assumption we have had to make, beyond our orderings, is that any stochastic component of the magnetic field (if present) satisfies δB stoch /B δ.

The Collisional Limit
The equations for δf e in Section 3.5 and Section 4.4 are derived under the assumption that ν ee ∼ ω ∼ k v th i . In many cases of interest, the collision frequency is in fact much larger than this: ν ee ω. We can take this limit as a subsidiary expansion, subsidiary to our expansion in the mass ratio. + Formally, we introduce the parameter ν * e = ν ee qR/v the , where q(ψ) is the safety factor (see Section 7.3 of [1]) and R the cylindrical radial coordinate (the major radius of the tokamak). As our initial ordering was ν ee ∼ ω, our ordering for ν * e is ν * e ∼ δ; for our collisional ordering we let ν * e δ and expand in The results of this expansion will be presented in the following sections. + Strictly speaking, to treat this as a subsidiary expansion we also require that ν ee k v the , so that it does not change the form of (43) and (46). However, if ν ee is large enough that ν ee ∼ k v the , then the homogenous equation for g e is w b · ∇g e = C[g e ]. The only solution of this equation is a perturbed Maxwellian with perturbed densities and temperatures that are functions ofψ only. Thus, the results of this section still hold.

Maxwellian Electrons
In our collisional expansion, C [g e ] is the dominant term in (92). Thus, up to corrections that are small in ω/ν ee 1, g e is given by Hence, g e is a perturbed Maxwellian to leading order: where δn H and δT e are arbitrary functions of r -δn H will shortly be shown to be the homogenous (i.e., a function only ofψ) part of the electron density perturbation and δT e the perturbed electron temperature. Now, we know from Section 3.1 that g e is a solution of (49). Substituting this into (49), we see that δn H /N e and δT e /T e are functions ofψ only (the derivatives in (49) being taken at constant ε e ). Thus, in the this limit, collisions trap and detrap electrons rapidly enough that there is no longer a distinction between passing and trapped electrons (α-dependence being the signature of kinetic trapped electrons) -all electrons sample the entire flux surface on the fluctuation timescale.
Let us now further simplify our solution for g e . Substituting (115) into (81), we obtain From this expression, we observe that we were justified in our notation for δT e -it really is the perturbed electron temperature relative to the moving surfaces. We also note that ζ and δn H both contribute to the perturbed density. As we have the freedom to add an arbitrary function ofψ to the definition of ζ -it does not change the velocity of the flux surface (see the end of Section 3.2), we redefine ζ to be ζ − T e δn H /eN e and so eliminate δn H entirely. Thus, our solution for δf e becomes where the entire density perturbation is contained in the Boltzmann-like first term.

Fluctuating Continuity Equation
We now derive an equation for ζ in the collisional limit.
Using the solution (115) for g e in (94), we obtain, to lowest order in δ/ν * e , ∂ ∂t ψ ,α,l where we have used the fact that, as we have absorbed the density perturbation into ζ, d 3 wg e = O((δ/ν * e ) n e ). As in Section 4.1, this is just the fluctuating continuity equation for electrons.

Fluctuating Energy Equation
To close our system of collisional equations, we now need the fluctuating energy equation from which we will determine δT e . Thus, we multiply (92) by m e w 2 /(2T e ) − 3/2 and integrate over all velocities to find where q e = d 3 w (m e w 2 /2T e − 3/2) h (1) e . As we do not know q e , we wish to eliminate it from this equation. Therefore, we average (119) over the perturbed flux surfaces and obtain ∂ ∂t where we have used (G.6) to eliminate some terms. Let us examine the terms in (120). Other than the last term on the right-hand side, which is due to the viscous heating of electrons, these terms seem fairly opaque. However, they all vanish if we apply the low-Mach-number ordering from Section 11.1 of [1] -we can move n e outside the flux-surface averages and then use the properties of the flux-surface average and the drift velocities to eliminate the terms in question. Thus, we conclude that these terms, which disappear if the Mach number is small, are to do with the exchange of potential energy of the electrons with thermal energy. Equation (120) completes our solution for δf e . In the next section, we close the system by working out the electron contributions to the field equations in the collisional limit.

Field Equations
We now use our solution for δf e to simplify the field equations. Substituting the solution (117) into (111) gives the following equation for δϕ : Similarly, by substituting (117) into (112), we obtain an equation for δB in the collisional limit: The field equations (121) and (122) close the collisional electron model, which we now summarise.

The Collisional Electron Model
The collisional electron model is a simplified version of the model given in Section 5. We have replaced the distribution functions g pe and g te by a single perturbed electron density and temperature. The complete set of equations describing this model is: • (67) and (68) to evolveψ andα which define the moving coordinate system, • ζ is given by (118), • δT e is given by (119), • (121) and (122) determine the fields δϕ and δB respectively, • and the ion dynamics (h i ) are given by (24).

The Transport Timescale
The procedure of Section 5 is sufficient to calculate the evolution of our plasma on the fast (turbulent) timescale. Our expansion in δ also has consequences for the evolution of n e , T e and ω on the slow (transport) timescale. In Section 8 of [1], transport equations for n e and T e were derived by considering the transport of particles and energy across ψ = const. surfaces. It is clear from the results of Section 3.2 and Section 3.3 that, for electron transport, it is more appropriate to considerψ = const. surfaces.
To this end, in Appendix H, we rewrite the electron kinetic equation relative to the moving flux surfaces. We then average this equation over the turbulence and over the fluctuating flux surfaces. Finally, we integrate over velocity space to find transport equations for particles and multiply by ε e and integrate over velocities to obtain the electron heat transport equation. This procedure is entirely analogous to that undertaken in Section 8 of [1] save for the extra complexity introduced by the rapid fluctuation of theψ = const. surfaces. We summarise the results of these calculations in the following sections.

Particle Transport
In Appendix H.3, we derive the electron particle transport equation, correct to lowest order in δ: where S (n) e is defined by (164) and the particle flux is given by Note that whilst we use the same notation for the particle flux Γ e as in [1], it is not the same quantity. In this paper the transport equations are derived for transport across ã ψ = const. surface not a ψ = const. surface. Several important features of this flux immediatly present themselves. Firstly, only the trapped electrons contribute -the passing particles play no part in the transport. Secondly, there is no "flutter" transport (transport due to correlations between δu and δB ⊥ despite the fact that the contribution to (170) of [1] from δA is clearly not zero. This is, in fact, the primary advantage of writing our transport equations with respect to the moving surfaces, we observe this cancellation analytically rather than having to recover it via a very careful evaluation of (166) and (170) of [1]. Finally, we see that drifts that do cause transport are the drift due to ζ and the ∇B drift due to δB . This supports our interpretation of the drift due to ζ being the drift of electrons across the moving surfaces.
7.1.1. Electron Transport in the Collisional Model As we proved in Section 115, in the collisional limit the electrons are trapped and de-trapped fast enough to remove the distinction between trapped and passing particles -in effect all electrons are passing. This immediatly implies that if the electrons are collisional, there is no transport of electrons. Thus, the electron transport time is O(δ −1 (ω/ν ee )) longer than the ion transport time -collisions are actually beneficial. Naturally, with sufficiently large collisions, neoclassical and classical electron transport reappear in the fluxes and the confinement worsens again.

Electron Heat Transport
Similarly, in Section Appendix I.4, we derive the electron heat transport equation, to lowest order in δ: and the heat sources are given by and P turb e = − d 3 w g e ∂ ∂t ψ ,α,l eζ − µ e δB turb ψ + e n e ζ∇ · u eff turb ψ .
In these expressions for the heating, the divergence of u eff is given by This expression is correct and closed to the required order in and δ. Again, it is crucial to note that whilst we use the same notation for q e , P comp s , P pot s , and P turb e as in [1], these are not the same quantities.
Let us now interpret the heat flux and the heating terms. The heat flux parallels the particle flux: there is no flutter transport, only trapped particles contribute to the flux, and they contribute via the E × B-drift across flux surfaces and the ∇B-drift due to δB .
Of the heating terms, the first, P comp s is just the compressional heating due to the mean motion of the flux surfaces. The second, P pot s , is the change in the thermal energy of the electrons due to the change in their electrostatic potential energy. * Finally, we have the turbulent heat source P turb e . We interpret the terms in the expression (129) for * For further discussion of the nature of this potential-energy-exchange term, see Section 8.3.5 of [1].
P turb e as follows: the first term is heating due to work done against the parallel electric field (through ζ); the term involving the time-derivative of δB is in fact composed of the electron viscous heating, compressional heating due to the fluctuating motion of theψ = const. surfaces and work done bending the magnetic field. This can be seen by substituting for δB from (130) -the term on the left-hand side of (130) is the compressional heating, the second term on the right is part of the electron viscous heating. The final term is the potential energy contribution to the compressional heating due to the fluctuating motion of the exact flux surfaces.
where ∇ · u eff is given by (130) as before. Thus, despite the lack of any transport, the electron response to ion-scale turbulence can still dissipate free energy. All other terms in (125) remain unchanged.

Momentum Transport
As the momentum transport involves all species, not merely electrons, we do not write a new equation for it in the moving frame. The transport equation for toroidal angular momentum is (Equation (179) where the moment of inertia of the plasma is now given by and the total momentum flux is where the sums over s = i are taken over all ion species and fluxes π (ψφ) s and Γ s for the ions are given by (184) and (170) of [1], respectively. The electron momentum flux is now (see Appendix H.5) The three evolution equations (123), (125) and (132) close our system on the transport timescale. The fluxes are written entirely in terms of fluctuating quantities that our reduced electron model provides (note that nowhere do we require the part of g pe that is odd in σ).

The Collisional Limit
In the collisional limit, g e is Maxwellian, and so isotropic in velocity space -see (115). Therefore, the velocity integral in (135) vanishes and so there is no angular momentum transport due to the electrons.

Summary
In this section we summarise our work, first by comparing it in detail to previous models of the electron response, then by briefly collecting our results and outlining future applications of this work.

Comparison to Previous Models
Many previous works have introduced physically-motivated model responses for the electrons. Some are based on the linear electron response [10] and others use fluid moments of the electron distribution function to try and capture the passing particle response [11,14,27]. Our model improves upon this by being a rigorous consequence of the electron gyrokinetic equation, clearly demonstrating how the electron response fits into the full turbulent transport framework. Also, by retaining the bounce-and fluxsurface-averaged kinetics of the electrons, we retain the distinction between the passing and trapped populations and retain drift resonances that may be lost in a fluid picture. We recover some of the results of the fluid models in the collisional limit, where there are no trapped electrons. In particular, the fact that δT H is only a function ofψ means that the temperature has equilibrated along perturbed field lines as expected. Even in this limit, the previous models did not allow for a rapid toroidal rotation (comparable to the ion thermal speed). Also all these models neglect δB and only retain δA in the perturbed magnetic field.
One of the original and still widely-used models for the electron response is the adiabatic electron response corrected for zonal perturbations [28]. It can be shown that this is a rigorous limit of our model in the limit of vanishing β (the electrostatic limit) and collisional electrons. For collisionless electrostatic turbulence, similar equations have been derived rigorously in [9]. Our model extends this to include fully electromagnetic turbulence and rapid toroidal rotation. The electrostatic model of [9] also neither handles general magnetic geometry nor discusses the behaviour of passing electrons. This limit is useful as an easily-implemented stepping stone between the adiabatic electron model (strictly only valid for quite collisional plasmas) and simulating the full electron gyrokinetic equation (which is either computationally prohibitive or formally invalid depending on whether one chooses to resolve the electron scales or not).
The most complete previous model is that of [12]. However, this model is derived under two restrictive assumptions. Firstly, the mass ratio is assumed to be much smaller than we assume here (δ ∼ rather than 1 δ which we assume) and secondly δB is assumed to be identically zero. We note that (81) corresponds directly to equation (54) of [12], but that the further assumption made in [12] of (in our notation) g pe = 0 is unjustified (g pe = 0 is not a solution of (110)) and hence neglects the electron response that is zonal with respect to the perturbed flux surfaces. In addition, the model in [12] does not include the effects of sonic toroidal rotation.
Finally, we must discuss how our model compares to the common practice of simulating the full electron gyrokinetic equation at ion scales (so-called "kinetic electrons"). Firstly, the electron gyrokinetic equation clearly contains electron-gyroscale and electron-bounce-frequency physics so simulating only the ion scales is formally invalid -deliberately leaving the simulation unresolved in some sense. Secondly, how an unresolved simulation of the electron gyrokinetic equation behaves will be sensitive to the numerical algorithm used -it is possible to design a numerical algorithm that numerically bounce-averages the kinetic equation in the long-timestep limit, but of the commonly-used gyrokinetic codes only GS2 [29] does this. Indeed, the only way to confirm that simulating "kinetic electrons" gives the correct physics is to either compare with full two-scale simulations that resolve the electron dynamics as well as the ion dynamics or to compare with simulations of the electron model presented in this paper. As well-resolved two-scale simulations are generally beyond current computational resources, simulations of our model are the only way to confirm the validity of this practice.

Conclusions
In this paper we have derived equations describing the electron response to ion-scale turbulence. The fast electron timescales have been eliminated from this system. The scale separation between electron and ion scales both motivates and enables the derivation of such a model. Starting from electron gyrokinetics [1], we systematically expand our equations in the square root of the electron-to-ion mass ratio to obtain our model. As discussed in the introduction, as studies of turbulence at non-zero β become more and more important -both for conventional and spherical tokamaks -having an electron response model that handles this rigorously will be necessary. We have derived full (β ∼ 1) collisional and arbitrary-collisionality models -summarised in Section 6 and Section 5 respectively.
The procedure by which we constructed this model, a rigorous expansions in m e /m i 1, allowed us to clearly delineate which effects are due to electron-scale physics and which to ion-scale. Most importantly, the elimination of the electron scales has a profound impact upon the structure of the fluctuating magnetic field. We prove in Section 3.2, that purely ion-scale turbulence conserves the total magnetic flux, and, moreover that a velocity field can be found into which the magnetic field is frozen.
This generalisation of the usual MHD frozen-in theorem demonstrates that it is only electron-scale effects which can change the magnetic topology. The inclusion of β ∼ 1, trapped particles, sonic rotation, and arbitrary axisymmetric geometry do not change this fundamental fact.
Similarly, whilst deriving transport equations across the exact flux surfaces, we have obtained another general property of ion-scale turbulence: passing electrons do not contribute to transport fluxes (see Section 7). This result could prove useful in designing magnetic geometries that have beneficial trapped-particle properties, safe in knowledge that the passing particles never contribute to the transport.
In addition to these general results, two other ancillary results have been proved in the course of this work. Firstly, in Appendix C, we proved that the topology of the mean magnetic field (i.e., the q(ψ) profile) evolves on the resistive rather than the transport timescale. Importantly, this result holds independently of the nature of the turbulence. Secondly, in Appendix J, we prove that the commonly-used adiabatic electron model is, in fact, the correct model for collisional (ν ω) and electrostatic (β 1) electrons. Throughout our derivation, we have simply assumed that the magnetic stochasticity produced by electron-scale fluctuations is sufficiently small that we can neglect it. This assumption must be investigated in detail and will be explored in future work. Similarly, the potential effects of small resonant regions around rational surfaces (which would give a finite width to eachψ-constant surface) should be examined closely. Indeed, a natural extension of this work would be to include the reaction of a small amount of purely electron-scale turbulence upon the ion scales. As always, the final arbiter of the usefulness of these models will be their numerical implementation and application to simulating moderate-and high-β ion-scale turbulence, in e.g. nonlinear studies of finite-β ITG turbulence or kinetic ballooning mode turbulence.

Acknowledgments
Part of this work was carried out at meetings supported by the Wolfgang Pauli Institute and the Leverhulme Trust Academic Network for Magnetized Plasma Turbulence. During this work, IGA was supported by a CASE EPSRC studentship jointly with the EURATOM/CCFE Fusion Association and by a Junior Research Fellowship at Merton College, Oxford. This work was carried out within the framework of the European Fusion Development Agreement, the views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A. Derivation of (42)
The gyrokinetic potential χ is Using (16), we perform a Taylor expansion about r = R e to find where we have used the fact that w · δA/c δϕ . Gyroaveraging this expression, the third and fourth terms vanish to leave where we have used w ⊥ w ⊥ R = w 2 ⊥ /2(I − b b), with I the unit dyadic, and δB b = ∇ × δA ⊥ .

Appendix B. Derivation of (57)
Expressing E in terms of potentials, we arrive at the following expression for E · b: Manipulating the first term of this expression, we find that where we have used (6), (9), and (12).
where in the last line we have used (45), u · ∇b = b · ∇u ((A.17) of [1]), and neglected δB · (∇u) · δA as it is O( The final term in this expression is O( 3 (v th i /c) B) and can be neglected. We can then use δϕ = δϕ − δA · u/c to finally write the parallel electric field as

Appendix C. Derivation of (58) and Conservation of Mean Magnetic Flux
In this Appendix, we first derive (58), as required for the derivation in Section 3.2. Then we prove that there is a field-line-preserving velocity for the mean magnetic field, irrespective of the fluctuations. Finally, we demonstrate that this implies that there are two timescales for the evolution of the mean magnetic field: the transport timescale τ E on which the flux ψ(R, z, t) evolves and the resistive timescale τ η ∼ δ −1 τ E on which the safety factor q(ψ, t) evolves. We now proceed to derive (58) in a way that exactly parallels the derivation of (54) from (24). The neoclassical drift-kinetic equation that determines F 1e is ((112) of [1]) The first term on the left-hand side and the second term on the right-hand side are O( 2 δ −1 Ω i F 0e ), all other terms are O( 2 Ω i F 0e ). Thus, to lowest order in δ, we obtain where we have dropped the distinction between R e and r as it is small in . This equation is completely similar to the first line of (46) and so we seek solutions in an identical fashion. Let (c.f. (48)) F 1e = λ(r, ε e , µ e , σ)F 0e + G e (r, ε e , µ e ) + O ( δF 0e ) , where G e satisfies w b · ∇G e = 0. (C.4) As G e is axisymmetric and so cannot depend on α, we immediatly solve (C.4) to find G e = G e (ψ, ε e , µ e ). Substituting (C.3) back into (C.2) and dividing by F 0e , we obtain By the same logic as in Section 3.1, λ is a function of r only (this time there is no gradient term to drive a temperature perturbation). Moreover, by absorbing any density moment of G e into λ, we can assume that d 3 wG e = 0 and so where n 1e is the first-order correction to the mean electron density: n 1e = d 3 wF 1e . Using this result in (C.5), we finally obtain This is precisely (58) as required in Section 3.2.
We now go to next order in δ to obtain an equation for G e . To this order, we have 1e is the next correction to F 1e . Using (116) of [1], we have Multiplying this equation by G e /F 0e , integrating over all velocities, and flux-surface averaging, we obtain The only solution to this equation is for G e to be a Maxwellian with a density and temperature that are functions of ψ only. We can absorb any such solution into F 0e . Thus, G e vanishes and the solution for F 1e is merely Substituting (C.7) into the expression for E · B in terms of the potentials ϕ and A, we have This form of the mean electric field is entirely analogous to (55). Thus, by the arguments of the second paragraph of Section 3.2, the mean magnetic flux is conserved and there is a field-line-preserving velocity for the mean magnetic field. Finally, we look at the implications of this for the evolution of the mean magnetic field given by (6). This field is determined by the two functions I(ψ, t) and ψ(R, z, t), which a priori we expect to evolve on the transport timescale τ E . In [1], it is shown that instead of I, we can use the safety factor q(ψ) given by to characterise the magnetic field. It is this quantity that will evolve slowly in the limit of δ 1. The evolution of q(ψ) is given by [1]: (C.14) Using the expression (C.12) for E · B, we see that and so the right hand side of (C.14) is O( 3 δΩ i ) = O(τ −1 η ), and thus q evolves on the resistive timescale, not the transport timescale. Since q(ψ) is the number of times a field line winds around the vertical symmetry axis for one poloidal transit around the magnetic axis, a topological quantity, it changes only when the magnetic topology changes. However, the shape of the flux surfaces changes on the transport timescale in response to changes in the pressure profile (via the Grad-Shafranov equation, see Section 7.3). It is therefore clear from (C.13) that I(ψ) also changes on the transport timescale.

Appendix D. Representations of δB
Once we introduce the Clebsch form of B we have two equivalent expressions for δB: firstly, the usual expression in terms of δA and δB , and secondly the expression in terms ofψ andα, Let us investigate this second expression for δB. Writing the gradients ofψ andα in terms of the (ψ,α,l) coordinates we have Examining the size of the various terms in this equation we see that becauseψ andα can vary on the short perpendicular spatial scale (∇ ⊥ ∼ ρ −1 i ). Thus, the first term in (D.3) is comparable to B. Similarly, we have that as parallel variation is much weaker than perpendicular. Thus, the second and third terms in (D.3) are smaller than the first by one power of . In order that δB is indeed smaller than B we require that Thus, to lowest order,ψ andα are not "independent". This represents the fact that, to lowest order, u eff is an incompressible flow and thus cannot generate any δB . Taking the scalar product of (D.2) with ∇ψ we obtain Substituting for δB via (D.1), we are able to relateψ to δA : Similarly, we can find an expression for δA in terms ofα: Hence, given any one of δA ,ψ orα we can determine the other two from the appropriate pair of equations from (D.6), (D.8) and (D.9). We do not use (D.3) to determine δB since for this we must knowψ andα to O( 2 ) or equivalently the O( 2 v th i ) compressible corrections to u eff -and we do not. Instead we use (32) to determine δB since it involves only lowest order quantities. From the above results we can also form two expressions for ∇δA in terms ofα andψ: and We will require these expressions when we come to eliminate δA from our equations.

Appendix E. Derivation of (79)
In the definition (56) of u eff we expand E in terms of potentials, keeping terms up to O( v th i ), to obtain where we have used the fact that (∂/∂t + u · ∇) δA ∼ Ω i δA to rewrite the time derivative of δA in the second term. Using (12) in the first term in (E.1), we find where we have also used (6) and (45). The second tern in (E.1) can be rewritten as follows: Substituting (E.2) and (E.3) back into (E.1), we obtain where we have used b · u = Iω(ψ, t)/B and (45). Letting U (r) = b · u and substituting for ξ from (60), we find where we have used the fact that (c/B) b × ∇ (T e n 1e /en e ) ∼ 2 v th i and can thus be neglected. Finally, pickingξ = ln N e (ψ), we arrive at This is precisely (79), in which ζ is defined by (80).

Appendix F. Properties of the Bounce Average
In this appendix we prove the properties of the bounce average that were stated in Section 4.2.
First, as the second portion of the integration only differs from the first in the sign of σ and the direction of integration, we have that, for any function a that is independent of σ, w a(l,ψ,α, ε e , µ e , t) = 0. (F.1) Secondly, combining (65) and (99), we discover that, for any single-valued function a, w b · ∇a(l,ψ,α, ε e , µ e , σ, t) = w ∂a ∂l where the contributions from the endpoints has vanished as at the bounce points a(σ = +1) = a(σ = −1) in order that a is continuous through w = 0.
We now prove that V De · ∇ψ = 0. Starting from the formula for V De , (89), we can show that, to leading order in δ, where the gradient is taken at constant energy ε e and magnetic moment µ e . Using this form for V De , we have that By using the expression of a divergence of a vector field A: we can write (F.4) as where we have used (F.2). Now, 1 Thus, as w is axisymmetric and so independent of α and thusα (to lowest order). We can now take the bounce average of this expression and use (F.2) to get the required result. Finally, we prove that the bounce average commutes with the time derivative in the field-aligned coordinates to the required order in . The integration overl is taken at constant t,ψ,α, ε e and µ e so clearly commutes with the time derivative in these variables. The fact that the denominator in (99) varies slowly in time, follows from the fact that it is in fact just the time taken for an electron with a given location, energy ε e and magnetic moment µ e to travel along a field line from one bounce-point, to the other, and back again. As the δB ⊥ ∼ B and k ∼ l, the contribution to the distance along the perturbed field line from the fluctuations is small. As w varies slowly in space, its value along the perturbed field line is almost the same as its value along the unperturbed field line. Thus, the bounce time along the perturbed field is a slowly varying function.

Appendix G. Properties of the Flux-Surface Average
As the flux surface average is defined in exactly the same way as the average over equilibrium flux surfaces defined in [1] (the fact that the integral is taken at constant ε e and µ e rather than v does not affect the derivation), we can use the results of Section 3.4 of [1] with a simple change of notation (ψ becomesψ). Thus, for any function A(r, v) we have where the integration is taken over the surface labelled byψ and Similarly, for any vector field A(r, v) we have where the derivative are all taken at constant ε e , µ e and ϑ.
We now prove the results that are stated without proof in Section 4.3. Using (G.3), we can show that Secondly, we prove (108). Working from (F.3), we have However, as, to lowest order in , B and |w | are independent ofα, we have where we have used the fact that λ is a fluctuating quantity to neglect ∂λ/∂l compared to ∂λ/∂α. Inserting this into (G.8) and using (G.5), we obtain (108) as promised.
Finally, as we can write V D in the form we can prove that and in the same manner as (G.7) and (108), by using the fact that T e = T e (ψ) + O( T e ) which commutes with the flux-surface average. As the magnetic perturbation is small, the fluctuations cannot change the area of a flux surface by more than an a fraction of order O( ). Thus V can change by at most V in one fluctuation turnover time, and so ∂V /∂t ωV . This means that the flux-surface average commutes with the time derivative at constantψ,α,l:

Appendix H. Transport Equations for Electrons
In this Appendix we derive the transport equations presented without proof in Section 7. We proceed from the kinetic equation written in r and v variables It will be convenient to introduce the peculiar velocity w = v − u eff , at which point we can write (H.1) as (H.2) To explicitly evaluate fluxes, we will need to use (H.2) to evaluate the gyrophasedependent piece of f e correct to second order in . Taking (H.2), dropping terms smaller than O( 2 Ω i f e ), and dividing by e B/cm e , we have

(H.3)
In this equation, the second-order parts of f e appear only on the left-hand side and so we will use this to eliminate the second-order parts of f e from our calculation. The structure of this appendix is as follows: first, in Appendix H.1 and Appendix H.2 we introduce some averages that will be needed later. In Appendix H.3, we derive the particle transport equation for electrons across the moving surfaces. In Appendix H.4, we derive the transport equation for electron heat. Finally, we evaluate the flux of toroidal angular momentum due to the electrons.
Detailed derivations of the intermediate results used in this appendix are provided in Appendix I.
Appendix H.1. The Turbulence Average at constantψ: In order to derive the transport equations, we will need some new averages. We define a new average over the fluctuations where the perpendicular spatial integrals are taken between surfaces of constantψ andα rather than ψ and α. Similarly, the time integral in this average is taken at constantψ,α andl. We denote this new average by · turb : where the spatial integral is taken over a small patch that has dimension λ (the intermediate length scale) in each direction -for reference compare with (15) of [1].
This average commutes with time derivatives at constantψ and with flux-surface averages overψ-constant surfaces.
As usual, the composition of this average and the average overψ-constant surfaces is an average over time composed with an average over an annular region betweenψ =ψ 0 andψ =ψ 0 + ∆. It is important to note that this average differs from the turbulence average in [1] only by terms that are small in . This will be crucial, as we wish to replace all instances of this average by their unprimed version in the final equations.
Appendix H.2. Flux Surface Averaging on the Transport Timescale: We will want to average the transport equations over the fluctuating flux surfaces. For this we need two identies. The first is just a specific case of (G.3). If A(r, t) is a function only of space and time then where the derivative on the right is taken at constant t only.
Secondly we need to interchange flux-surface averages and time derivatives. We do this via which is proven in an identical fashion to (45) of [1] except that we retain the velocity explicitly rather than replacing it with ∂ψ/∂t via (67).

Appendix H.3. Particle Transport
We now move to the meat of the problem. Deriving our first transport equation.
Integrating ∼ n e /τ E . The particle flux Γ e is defined by (H.10) In Appendix I.3, we use the kinetic equation to find the following explicit expression for the particle flux: Appendix H.4. Energy Transport and Heating Multiplying (H.2) by m s w 2 /2, integrating over all velocities and averaging over fluxsurfaces and the turbulence we obtain 1 V where the electron thermal energy flux is defined by where V ψ is the velocity of the mean flux surfaces defined in (44) of [1], and the exchange between thermal energy and electrostatic potential energy: (H. 17) In (H.17) and (H.15), the divergence of u eff is given by This expression is correct and closed to the required order in and δ. In Appendix I.5, we show that the heat flux can be written in terms of g e as We now collate results from the rest of the paper to write down a convenient form for f e . First, using the general solution (13) for f e , we can write , ε e , µ e , ϑ, t) 1e ∼ δ f e contains the higher-order (in δ) contributions to F 1e . We also wish to write F 0e in terms of the new velocity variable w. To this end, we introduce the following quantity (I.4) Using (60) and (80), we obtain the following expression for ξ eξ T e = ln N e (ψ) + eϕ 0 T e + n 1e n e + eζ T e , (I. 5) in which T e is evaluated atψ and we have let ξ = 0. Substituting this into (I.4), and using the definition (21) of F 0e , we can show that 1e by some small gyrophase-independent pieces -we will never need F (H) 1e explicitly so can absorb any gyrophase-independent function into it. Because of the convenience of the form of F e , we choose to write f e as Using our solution (81) for δf e and the definition of h e , we have that, to lowest order in δ, We will use this equation later to eliminate h e from our final results.
Appendix I.2. An equation for f 2e As mentioned above, we will need properties of f 2e to calculate the fluxes in our transport equations. We now derive such an equation by using the solution (I.7) in (H.3) -we will need to retain corrections to f 2e up to O(δ 2 f e ). We perform this substitution stepby-step. The preliminary step is to notice that all the terms in the second line of (H.3) are O(δ 2 2 f e ) or smaller and so do not contribute in this calculation. Next, we compute the left-hand-side of (H.3) for F e : where we have used the definition of R e and the fact that b · ∇ψ = 0. Using this result and the definition (I.4) of F e , we can show that We can rewrite this in terms of V De as Substituting this back into (I. 15) and thence back into (I.12) we arrive at: +O( 2 δ 2 f e ).
(I. 18) Estimating f 2e from this equation, we see that f 2e ∼ δ 2 f e , rather than f 2e ∼ 2 f e which we would naïvely expect -this is because it is only the gyrophase-dependent piece of f 2e we need and, to lowest order in δ, f 2e is gyrophase-independent.
Using the solution (I. 7) for f e in the definition of the particle flux, we havẽ Substituting for F e from (I.4) and performing the velocity integration, we havẽ (I.20) Using (79) for u eff and using the definition of the flux-surface average in the first term, we haveṼ Integrating by parts with respect toα in the first term, we find that The calculation of the term involving h e is particularly involved. In Appendix I.6, we prove the following identity: Ṽ h e w · ∇ψ where we have integrated by parts in the term involving f 2e .
(I. 25) where we have integrated by parts with respect to w in the second line and used the gyrophase-independence of g e in the third line. We have also used the isotropy of F 0e and the gyrotropy of g e to eliminate terms in (I.18) that contain an even power of w ⊥ (we are multiplying by one power of w ⊥ so this becomes an odd power and thus vanishes upon integration over ϑ).
We now prove that the last line of (I.25) is small and can be dropped. Writing it in terms of µ e and using (F. 8 where we have been able to write (F.8) in terms of w because the difference between w and w is O( v th i ) and thus negligible. Finally, integrating by parts with respect toα in the second line of (I. 25) and dropping all the references to w, we have where the passing electron response has vanished because ∂g pe /∂α = 0 (at constant ε e and µ e ). This is (H.11) as required.
Appendix I.4. Derivation of the Heat Transport Equation As in Appendix C.4 of [1], it is convenient to adjust the form of the pressure evolution equation before calculating the heat flux. In particular, we subtract where we have integrated by parts with respect to time in the second term. The first of these terms cancels with the corresponding term in (I.29) and the second and third will become part of the turbulent heating. We will need to process terms involving ∇ · u eff several times in this section. In Appendix I.4.2, we show that we can write it in a closed (if unwieldly) form -see (I.45). For conciseness we will not expand ∇ · u eff in this section. Now, the time-derivative term is small: Estimating the size of the contribution from the gyrophase-dependent piece of f 2e ∼ δ 2 f e , we see that it is O(δ 3 Ω i n e T e ) and so we can neglect it. Now, as ∇ · u eff turb ∼ b · ∇u eff · b turb ∼ O( 3 Ω i ), we can also neglect the contributions from where we have used the fact that, to lowest order in , the flux-surface average and the turbulence average commute. As u eff − u ∼ v th i we can use (79) to express u eff in the second term. Doing this, we find that the entire second term is O(δ 3 n e T e Ω i ) or smaller and can thus be neglected. In the first term, we cannot use (79), but instead, we interchange the divergence and the turbulence average to find that − m e d 3 wv · ∇u eff · wF e ψ turb = −T e n e ∇ · V ψ ψ , (I.37) where we have used u eff turb = V ψ -the average of the velocity of the exact surfaces is just the velocity of the average surfaces. We now work out the corresponding term involving g e . Using the gyrophase-independence of g e , we have where the passing electron response has vanished because ∂gpe ∂α = 0 (at constant ε e and µ e ). This is (H.19) as required.
Appendix I.6. Derivation of (I. 23) In order to derive this result, we will need a concise notation for the following operation: This is an important operation as ∇l is not parallel to the exact field, but the strong perpendicular variation is contained inψ andα. We see that this term is just the ∇B drift in the exact magnetic field. Now we handle the first term on the right-hand side of (I.60). Integrating by parts with respect to space, we have: where we have used the fact that ∇ψ = ∇ ⊥ψ and the fact that we can then choose to take that derivative at constant ε e and µ e rather than at constant w. The crucial step in this derivation is to now writeψ inside the first bracket as ψ(R e ) + (r − R e ) · ∇ ⊥ψ . Doing this, we obtain and that typical fluctuating velocities scale with the ion thermal speed (rather than the Alfvén speed): and that this scaling also applies to the fluctuating E × B velocity: In order to consider electrostatic turbulence, we assume that δB ⊥ and δA are small: These orderings imply thatψ andα differ from ψ and α by amounts that are small in √ β.
Estimating the size of terms in (112) we see that δB B ∼ β , (J.5) and so we can drop all terms containing δB .
As we wish to obtain the simplest equations possible, we will also neglect all finite-Mach-number effects. Thus, we combine this expansion with the low-Mach-number expansion of Section 11 of [1]. The primary result of the low-Mach expansion is that we can neglect ϕ 0 and the poloidal variation of the mean density. In addition to this simplification, we can drop the distinction between δϕ and δϕ -from now on only working with the electrostatic potential δϕ.
In the subsequent sections we will apply this ordering to our collisionless and collisional equations. One result is common to both the collisional and collisionless electrostatic limits -the general solution for ζ. We obtain this by applying our electrostatic ordering to (82). The left-hand-side of (82) is √ β smaller than the right; thus, to lowest order, we have b · ∇ (ζ − δϕ) = 0, (J. 6) where we have also used the fact that b = b + O( √ β). Solving for ζ, we obtain ζ = δϕ +ζ(ψ), (J.7) withζ an arbitrary function to be determined. Using this result in the expression for the field-line velocity u eff , we see that u eff is directed within the flux surfaces and merely moves one field line onto another without distorting them. Thus, this solution for ζ is consistent with our assumption that δB is small.