Radially local approximation of the drift kinetic equation

A novel radially local approximation of the drift kinetic equation is presented. The new drift kinetic equation that includes both ${\bf E} \times {\bf B}$ and tangential magnetic drift terms is written in the conservative form and it has favorable properties for numerical simulation that any additional terms for particle and energy sources are unnecessary for obtaining stationary solutions under the radially local approximation. These solutions satisfy the intrinsic ambipolarity condition for neoclassical particle fluxes in the presence of quasisymmetry of the magnetic field strength. Also, another radially local drift kinetic equation is presented, from which the positive definiteness of entropy production due to neoclassical transport and Onsager symmetry of neoclassical transport coefficients are derived while it sacrifices the ambipolarity condition for neoclassical particle fluxes in axisymmetric and quasi-symmetric systems.


I. INTRODUCTION
Effects of neoclassical transport 1-3 on plasma confinement are more significant in stellarator and heliotron plasmas than in tokamak plasmas because, in the former, radial drift motions of trapped particles in helical ripples enhance particle and heat transport due to nonaxisymmetry of the magnetic configuration. [4][5][6] Conventional calculations of neoclassical transport fluxes are done applying radially local approximation to solving the drift kinetic equation, in which v d ·∇f are often neglected as a small term of higher order in the normalized gyroradius parameter δ ∼ ρ/L. (Here, v d , f , ρ, and L represent the guiding center drift velocity, the deviation of the guiding center distribution function from the local Maxwellian equilibrium distribution, the gyroradius, and the equilibrium scale length, respectively.) However, in stellarator and heliotron plasmas, this v d · ∇f term is known to be influential on the resultant neoclassical transport because it significantly changes orbits of particles trapped in helical ripples. Therefore, at least, the E × B drift part v E · ∇f in v d · ∇f has been kept in most studies of neoclassical transport in helical systems. [7][8][9][10][11][12][13][14][15] Recently, it was shown by Matsuoka et al. 13 that the neoclassical transport is significantly influenced by retaining the magnetic drift tangential to flux surfaces in v d · ∇f for the magnetic configuration of LHD especially when the radial electric field is weak. However, as pointed by Landreman et al., 14 stationary solutions of the drift kinetic equation with radially local approximation used require additional artificial sources (or sinks) of particles and energy when the above-mentioned drift terms are retained. In this paper, a novel radially local drift kinetic equation, which includes both E × B and tangential magnetic drift motions, is presented. The radially local guiding center motion equations do not satisfy the conservation law of the phase-space volume while the full guiding center motion equations do. This fact causes the difficulty in obtaining the stationary solution of the local drift kinetic equation. However, the new local drift kinetic equation, which is written in the conservative form, has favorable properties for numerical simulation such that any additional terms for particle and energy sources are unnecessary for obtaining stationary solutions. In addition, it satisfies the intrinsic ambipolarity condition for neoclassical particle fluxes in axisymmetric systems as well as in quasi-symmetric helical systems. 16,17 The present work also treats interesting issues regarding the entropy production rate and Onsager symmetry 18,19 for neoclassical transport equations resulting from the new local drift kinetic model. The rest of this paper is organized as follows. In Sec. II, we consider the full drift kinetic model based on Littlejohn's guiding-center equations 20 without radially local approximation. Particle, energy, and parallel momentum balance equations are derived from the full drift kinetic equation. These balance equations are flux-surface averaged to confirm that they contain the second-order terms in δ, which represent neoclassical transport across flux surfaces. Also, expanding the distribution function about the local Maxwellian, we rewrite the drift kinetic equation to explicitly show that the thermodynamic forces defined by the background density and temperature gradients and the parallel electric field cause the deviation f from the local Maxwellian. In Sec. II, a new drift kinetic model is constructed by applying radially local approximation to Littlejohn's guiding-center equations with keeping E × B and tangential magnetic drift velocities. The new local drift kinetic equation for f is shown to be compatible with the stationary solution and to give intrinsic ambipolar particle fluxes for axisymmetric and quasi-symmetric systems. In Sec. IV, we present another radially local drift kinetic equation, from which the positive definiteness of entropy production due to neoclassical transport and Onsager symmetry of neoclassical transport coefficients are derived although this lo-cal drift kinetic equation no longer guarantees rigorously the intrinsic ambipolarity of neoclassical particle fluxes for axisymmetric and quasi-symmetric systems. Finally, conclusions are given in Sec. V.

II. FULL DRIFT KINETIC MODEL
A. Drift kinetic model based on Littlejohn's guiding-center equations We denote the guiding-center variables by (X, U, ξ, µ), where X represents the position vector of the guiding center, U the parallel velocity, ξ the gyrophase defined by the azimuthal angle of the gyroradius vector around the magnetic field line, and µ the magnetic moment. The Lagrangian for the guiding-center motion is given by Littlejohn 20 as where the Hamiltonian H is given by Here, Φ denotes the electrostatic potential. Using Eqs. (1) and (2), the guiding-center motion equations are derived as where Ω = eB/(mc), , and A * ≡ A + (mc/e)U b are used, and the guiding-center drift velocity V gc is defined by the righthand side of the equation for dX/dt. On the right-hand side of the equation for dU/dt in Eq. (3), the last term U b · ∇b · V gc is smaller than other terms by the order of δ = ρ/L where ρ and L represent the gyroradius and the gradient scale length given by L ∼ B/|∇B| ∼ Φ/|∇Φ|. The Jacobian for the guiding-center variables is written as where x and v denote the particle position vector and the velocity vector, respectively. Then, the conservation of the phase-space volume which can be proved by using Eqs. (3) and (4). The drift kinetic equation for the distribution function F (X, U, µ, t) is given by where the total time derivative is denoted by˙= d/dt. In the right-hand side of Eq. (6), C(F ) is the collision term and the additional term S is given to represent external particle, momentum, and/or energy sources if any.
Here, S is considered to be of the second order in δ. We can also treat effects of turbulent fluctuations by Eq. (6) if we regard the second-order additional term S as the ensemble average of the product of fluctuation parts in the electromagnetic fields and the distribution function as shown in Refs. 21 and 22 where the notation D is used instead of S to represent the term including the effects of turbulent fluctuations. Using Eq. (5), the drift kinetic equation can be rewritten in the conservative form as B. Particle, energy, and parallel momentum balance equations Multiplying Eq. (7) with an arbitrary function A(t, X, U, µ) which is independent of the gyrophase ξ and taking its velocity-space integral, the balance equation for the density variable d 3 v F A in the X-space is derived as and the velocity-space integral is denoted by d 3 v = 2π dU dµ D for gyrophase-independent integrands. For the case of A = 1, Eq. (8) reduces to the timeevolution equation for the density d 3 v F , In deriving Eq. (10), the conservation law, d 3 v C(F ) = 0, is used. However, it is noted that, if we use the collision operator obtained by the transformation from the particle coordinates to the guiding-center coordinates with finite-gyroradius effects taken into account, the velocityspace integral d 3 v C(F ) does not vanish but it becomes the opposite sign of the divergence of the classical particle flux as shown in Refs. 23-25. Here and hereafter, we assume that the expression of C(F ) is the same as that of the Landau collision operator given in the particle coordinates for simplicity so that d 3 v C(F ) = 0 is satisfied and the classical transport is neglected. We next consider the energy E = H [see Eq. (2)] as A in Eq. (8) and obtain the energy balance equation, where the total time derivative of the energy is written asĖ We easily see from Eq. (12) thatĖ = 0 for the stationary electromagnetic field. When we use the kinetic energy, another form of the energy balance equation is given by where the total time derivative of the kinetic energy is written asẆ The parallel momentum balance equation is derived from Eq. (8) with A = mU as We now use and rewrite Eq. (16) as Here, the density n, the parallel flow velocity u , the pressure tensor P, and the parallel friction force F are defined by whereẊ ⊥ ≡Ẋ − (Ẋ · b)b. Note that the pressure tensor P consists of the Chew-Goldbeger-Low (CGL) tensor 26 P CGL and the viscosity tensor π 2 of the second order in δ, where π 2 satisfies π 2 : I = π 2 : bb = 0 and the deviation of F from the local Maxwellian distribution is considered to be of O(δ).
It is well known that, if we use the original Boltzmann kinetic equation instead of the drift kinetic equation in Eq. (7), we can derive the momentum balance equation, where the Boltzmann kinetic equation is assumed to also contain the source term S. In Eq. (20), the particle flow nu, the pressure tensor P, and the friction force F are defined by nu = d 3 v F , P = d 3 v F mvv, and F = d 3 v C(F )mv, where, exactly speaking, F = F (x, v, t) represents the particle distribution function given by the solution of the Boltzmann kinetic equation and it has a gyrophase dependence that is not included in the solution of the drift kinetic equation. Comparing Eqs. (18) and (20), we see that Eq. (18) coincides with the parallel component of the exact momentum balance equation in Eq. (20) except that the former contains nmu·∂b/∂t and the non-CGL viscosity tensor expressed differently from the one in the latter.
We now consider general toroidal configurations, for which the magnetic field is written in terms of the flux coordinates (s, θ, ζ) as where θ and ζ represent the poloidal and toroidal angles, respectively, and s is an arbitrary label of a flux surface. The poloidal and toroidal fluxes within a flux surface labeled by s are given by 2πψ(s) and 2πχ(s), respectively. The derivative with respect to s is denoted by ′ = d/ds so that ψ ′ = dψ/ds and χ ′ = dχ/ds. Taking the fluxsurface average of the covariant toroidal component of Eq. (20) and making the summation over species, we obtain the expression for the radial current as 6 where the superscript s and the subscript ζ represent the covariant radial component and contravariant toroidal component given by taking the inner products with ∇s and ∂x/∂ζ, respectively, and the subscript a is used to explicitly show the particle species. Using the symmetry property of the pressure tensor P, we can show that, for axisymmetric toroidal systems, where P s ζ = ∇s · P · ∂x/∂ζ. In axisymmetric and quasiaxisymmetric toroidal systems, 17 we have Then, using Eqs. (22)-(24) and P = P CGL + π 2 , we find that, even for axisymmetric toroidal systems in the stationary state (∂/∂t = 0) with S = 0, the surfaceaveraged radial current does not vanish exactly due to the second-order viscosity tensor π 2 as shown by However, it is shown in Ref. 17 that (π a2 ) s ζ is a small quantity of O(δ 3 ) in axisymmetric systems with up-down symmetry (as well as in quasi-axisymmetric systems with stellarator symmetry) where all terms in the toroidal momentum balance equation given from Eq. (22) vanish up to O(δ 2 ). The same argument as above can be done for other quasi-symmetric systems such as quasipoloidally-symmetric and quasi-helically-symmetric systems if stellarator symmetry holds. On the other hand, in axisymmetric systems without up-down symmetry, is not guaranteed. Then, the ambipolarity condition a e a n a u s a = 0 is not automatically satisfied on the second order in δ because of the thirdorder radial particle fluxes (c/e a χ ′ V ′ )∂[V ′ (π a2 ) s ζ ]/∂s driven by the second-order shear viscosity tensor components (π a2 ) s ζ [here, it is useful to formally regard the electric charge as the O(δ −1 ) quantity 27 so that the radial current due to the third-order radial particle flux is immediately found to be of the second order]. However, even in this axisymmetric but up-down asymmetric case, the second-order radial neoclassical particle fluxes driven by the CGL tensors still automatically satisfy the ambipolarity condition for the radial current up to the first order. 1-3

C. Drift kinetic equation expressed in terms of flux coordinates
Using the flux coordinates (s, θ, ζ), the drift kinetic equation, Eq. (6), is rewritten as In Eq. (27), the functions s(X, t), θ(X, t), and ζ(X, t) are defined by the inverse of X = X(s, θ, ζ, t), where t is generally included as a parameter. Denoting the Jacobian for the flux coordinates (s, θ, ζ) by the conservation law of the phase-space volume, Eq. (5), and the conservative form of the drift kinetic equation, Eq. (7), are rewritten as and respectively. For an arbitrary function A(s, θ, ζ, U, µ, t) which is independent of the gyrophase ξ, the phase-space integral is written as where represents the flux-surface average and denotes the radial derivative of the volume V (s) enclosed within a flux surface labeled by s. We now integrate Eq. (30) with respect to the coordinates (θ, ζ, U, µ) to obtain wherė The time-evolution equation for the surface-averaged density For the cases of A = W , Eq. (34) reduces to the surfaceaveraged energy balance equation, whereẆ is given by Eq. (15). In Eqs. (36) and (37)

D. Expansion about a local Maxwellian distribution
The zeroth-order solution F 0 of the drift kinetic equation, Eq. (26), is given by the local Maxwellian, which annihilates the collision term, The total time derivative of F 0 is written as where the zeroth-order density n 0 and temperature T 0 are flux-surface functions independent of (θ, ζ), and their time dependence follows the transport ordering, ∂/∂t = O(δ 2 ). The parallel electric field E is given by where Φ = Φ − Φ . We now define the first-order distribution f by where l dl represents the integral along the magnetic field line. Then, substituting Eqs. (41) and (42) into Eq. (26) yields where the thermodynamic forces are defined by  45) is of the second order in δ and this gives rise to global or finite-orbit-width effects on neoclassical transport.

III. RADIALLY LOCAL APPROXIMATION
Under the radially local approximation made here, the guiding center equations are written as where the second-order part −∇ Φ−c −1 ∂A/∂t of the electric field E is neglected and the guiding center drift velocity in Eq. (3)  In Eq. (46), the magnetic moment µ is allowed to vary in time such that conservation of the kinetic energy of the particle W = mU 2 /2 + µB, is satisfied. It might appear that the energy E = W + eΦ should be conserved instead of W . However, using Φ ≃ Φ , we find that the difference eΦ between E and W is approximately constant along the radially local guiding center orbit and accordingly the conservation of W is reasonable under the radially local approximation. We now define (V (rl) gc ) ⊥ by removing the radial component from (V gc ) ⊥ as where B * and E * in the definition of V gc given by Eq. (3) are replaced with their lowest-order parts B and −(dΦ/ds)∇s, respectively, and the factor α(Λ) is introduced to satisfy the condition (V (rl) gc ) ⊥ (µ = 0) = 0. Here, the ratio of the magnetic moment µ to the kinetic energy W = mU 2 /2 + µB is used to define the dimensionless parameter, Λ ≡ µB max /W , where B max is the maximum value of B on the flux surface. This parameter Λ is a measure for classifying the guiding center motion into either passing or trapped orbit. As Λ increases from 0 and approaches to 1, the orbit changes from the passing to the trapped one. Then, we assume that while α(Λ) = 1 except for an interval, 0 ≤ Λ < Λ 0 , where Λ 0 (≪ 1) is a small positive constant value. For example, α(Λ) is defined by We should note that influences of the magnetic and E × B drift motions are significant mainly for precession drift orbits of trapped particles, and that particles in the region, Λ < Λ 0 , are passing ones whose orbits almost coincide with field lines. Therefore, even if the functional form of α(Λ) and the value of Λ 0 are changed, the artificial reduction factor α(Λ) for Λ < Λ 0 is expected to cause little change in resultant passing particles' orbits except that the limiting condition, lim Λ→+0 (V (rl) gc ) ⊥ = 0, is rigorously satisfied. However, this insensitivity to the form of α(Λ) remains a future subject to be verified by numerical simulations.
It also should be mentioned that the radially local approximation described by Eqs. (46) and (48) is independent of what poloidal and toroidal angles are chosen for the flux coordinates. This is a favorable property that is lost in Ref. 13. We see that the radially local guiding center equations given by Eqs. (46) and the Jacobian D = B * /m for the phase-space coordinates (X, U, ξ, µ) [see Eq. (4)] do not satisfy the conservation law of the phase-space volume as shown in Eq. (5). This violation of the phase-space-volume conservation occurs even if B * is used instead of B in the denominator on the righthand side of Eq. (48) to define (V (rl) gc ) ⊥ . In the next section, we consider another Jacobian in order to recover the conservation law although, in this section, a simpler approximate Jacobian D 0 ≡ B/m is used. Also, we hereafter employ (X, W, U, ξ) as phase-space coordinates. Then, from the Jacobian D 0 ≡ B/m for (X, U, ξ, µ) with µ = (W − 1 2 mU 2 )/B, the Jacobian for (X, W, U, ξ) is derived as 1/m, which is constant in the phase space.
Using V (rl) gc , dU/dt, and dW/dt = 0 given by Eqs. (46) with (48) under the radially local approximation, the drift kinetic equation for the first-order distribution function f (X, W, U ) in the stationary state is written as The radial component of the guiding center drift velocity V s gc on the right-hand side of Eq. (51) is given by In deriving Eq. (52) from the guiding center drift velocity given in Eq. The fact that the Jacobian is constant is used in deriving Eq. (51) which is rewritten by using the flux surface coordinates (s, θ, ζ) as Here, we should note that, in Eq. (53), partial derivatives of the first-order distribution function f are taken only with respect to the three variables (θ, ζ, U ) and that the radial coordinate s and the kinetic energy W enter f (s, θ, ζ, W, U ) as constant parameters.
Taking the velocity-space integral of Eq. (51) yields the continuity equation in the stationary state, The parallel and perpendicular particle fluxes in Eq. (54) are defined by where the velocity-space integral is written in terms of the variables W and U by The diamagnetic flow Γ ⊥1 = n 0 u ⊥1 and the parallel flow Γ are of the first order in δ = ρ/L while Γ (rl) is of the second order. The collisional particle conservation law, d 3 v C L (f ) = 0, is used to obtain Eq. (54). Also, it should be noted that the boundary condition, (V from Eq. (51). We find that the flux surface average of the left-hand side of Eq. (54) automatically vanishes so that no particle source is required for obtaining the stationary solution. Thus, the radially local approximation presented here has self-consistency with neglecting the radial transport that causes variation in the surfaceaveraged particles' number [see Eq. (36)].
Next, we multiply Eq. (51) with (W − 5T /2) and take its velocity-space integral to derive where the parallel and perpendicular heat fluxes are given by and the collisional heat generation is defined by In Eq. (57), q which represents the collisional heat exchange balance that needs to be satisfied in the stationary state. Unequal temperatures T a0 = T b0 can occur in the case of m a /m b ≪ 1 or m a /m b ≫ 1 where the characteristic time of the collisional thermal equilibration between the species a and b is much longer than the 90 • scattering times due to like-species collisions characterized by τ aa and τ bb . Then, C ab (f a0 , f b0 ) does not vanish even for the local Maxwellian distribution functions f a0 and f b0 given by Eq. (38) and it describes the above-mentioned slow collisional thermal equilibration although the linearized collision operator C L used for the Eq. (53) does not include this equilibrium part of the collision term. However, the heat generation Q ab , which is defined by Eq. (59) with the linearized operator C L ab for collisions between different species a and b, generally remains nonzero (even for the case of T a0 = T b0 ). Therefore, Eq. (60), which is rewritten as Q a ≡ b =a Q ab = 0 (recall Q aa ≡ 0), is considered to be the physically reasonable condition that should be satisfied in the multi-species stationary state of the radially local model without requiring additional heat source or sink.
Multiplying Eq. (51) with mU and taking its velocityspace integral give the parallel momentum balance equation, where the first-order pressure p 1 and the viscosity tensors π 1 and π (rl) 2 are defined by and the parallel friction force is given by The first-order viscosity tensor π 1 is written in the form of the traceless part of the CGL pressure tensor as π 1 = (p −p ⊥ )(bb− 1 3 I), where p and p ⊥ represent the parallel and perpendicular pressures, respectively. On the other hand, the second-order viscosity tensor π (rl) 2 , which is given by the correlation between the parallel velocity U and the perpendicular drift velocity (V (rl) gc ) ⊥ , can not be written in the CGL form. We now multiply Eq. (61) with the magnetic-field strength B and take its magneticsurface average to derive which is used later to derive an alternative expression for the neoclassical particle flux. The radial neoclassical particle flux is written as where V s gc = V gc · ∇s is given by Eq. (52). Derivation of Eq. (65) uses the following formula, and the Boozer coordinates (s, θ, ζ), 28 for which the contravariant poloidal and toroidal components, B θ and B ζ , of the magnetic field B are flux-surface functions. We see from Eq. (65) that the neoclassical particle flux is caused by the spatial gradients of the first-order pressure and viscosity tensor. It can be shown that the second-order viscosity tensor π (rl) 2 defined in Eq. (62) satisfies In deriving Eqs. (67) and (68), it is convenient to write π (rl) 2 = A(Bw + wB) with w ≡ (b × ∇s)/B. Then, we find w · (∇ · π (rl) 2 ) = w 2 B · ∇A + A(B · ∇w · w − w · ∇w · B) = w 2 B · ∇A − A∇s · (∇ × w) = ∇ · (A∇s × w) and (b/B)·(∇·π It is found from Eq. (65) and Eqs. (67)-(69) that the second-order viscosity tensor π (rl) 2 in the radially local approximation cannot contribute to the neoclassical transport like the first-order pressure p 1 and viscosity tensor π 1 .
We now use Eqs. (61) and (68) to rewrite the expression of the radial neoclassical particle flux in Eq. (65) as Here, the first surface-averaged part on the right-hand of Eq. (70) represents the nonaxisymmetric part of the neoclassical particle flux, 6,19 Then, we find from Eqs. (70) and (71) that the radial electric current is written as where the quasineutrality a n a0 e a = 0 and the collisional momentum conservation a F a = 0 are used. For axisymmetric and quasi-axisymmetric systems, we have from which and are derived. Here, the quasi-axisymmetry means that the magnetic field strength B = |B| is independent of the toroidal angle ζ.
automatically in axisymmetric and quasi-axisymmetric systems. The intrinsic ambipolarity can be proved in the same way for all other quasi-symmetric systems such as quasi-poloidally-symmetric and quasi-helicallysymmetric systems. It is seen from Eqs. (22)-(25) that the radial current is closely related to the toroidal viscosity or the radial transport of the toroidal momentum. As remarked after Eq. (25), the ambipolarity condition is not guaranteed on the second order in δ for axisymmetric systems without up-down symmetry (as well as quasi-axisymmetric systems without stellarator symmetry) because of the component (π a2 ) s ζ of the second-order non-CGL viscosity tensor. We also note that the radial neoclassical particle flux defined by Eq. (65) is the second-order flux driven by the first-order CGL tensor which becomes a dominant part for nonaxisymmetric systems although it does not contain the third-order flux due to the second-order tensor. Therefore, Eq. (76) should be interpreted to imply that the intrinsic ambipolarity condition for the axisymmetric case can be correctly treated only up to the first order by the present radially local approximation. On the other hand, the second-order neoclassical radial flux (π a2 ) s ζ of the toroidal momentum in the axisymmetric but up-down asymmetric case can also be evaluated using the solution f of the radially local drift kinetic equation even without resort to the radially global model. This can be done by substituting the solution f into the formula for the toroidal momentum transport flux given by Eq. (18)

IV. ENTROPY PRODUCTION RATE AND ONSAGER SYMMETRY ASSOCIATED WITH NEOCLASSICAL TRANSPORT EQUATIONS
The neoclassical radial particle flux Γ ncl , heat flux q ncl , and parallel electric current J E = BJ / B 2 1/2 are defined in terms of the solution f of Eq. (51) or (53) by where the subscript a denotes the particle species. The linearized collision operator in Eq. (51) for the species a is defined in terms of the bilinear operator C ab for collisions between the species a and b by Here, C T ab (f a ) ≡ C ab (f a , F b0 ) and C F ab (f b ) ≡ C ab (F a0 , f b ) are called test-and field-particle collision operators, respectively, and they satisfy the adjointness relations, 30,31 (79) and Boltzmann's H-theorem, 30,31 Strictly speaking, the adjointness relations and the Htheorem are rigorously satisfied by the linearized Landau collision operator only for the case of T a = T b although they are still approximately valid even for 30 Since the drift kinetic equations for different particle species are coupled with each other due to the field particle collision operators, f a depends not only on thermodynamic forces (X a1 , X a2 , X E ) but also on those for b = a, (X b1 , X b2 ). Accordingly, we find that Γ ncl a , q ncl a , and J E in Eq. (77) are related to the thermodynamic forces through the neoclassical transport equations which are written as Here, the neoclassical transport coefficients (L 11 ab , L 12 ab , · · · ) are regarded as functions of the variables [E s (≡ −dΦ/ds), ∇s · ∇B, ∇s · (b · ∇b)] which characterize the perpendicular guiding center velocity (V (rl) gc ) ⊥ defined in Eq. (48).
We here examine the Onsager symmetry of the neoclassical transport coefficients. In order to prove the Onsager symmetry, the adjointness relations written in Eq. (79) and the phase-space-volume conservation along the collisionless guiding center orbit are required as shown in Refs. 19 and 29. However, in the radially local model based on Eq. (51), the latter condition ∇·V (rl) gc +∂(dU/dt)/∂U = 0 is broken so that the Onsager symmetry is not satisfied.
As noted before Eq. (51), the Jacobian for the phasespace coordinates (X, W, U, ξ) is given by 1/m. Here, we consider a modified Jacobian, which differs from the one mentioned above by the correction term d * of O(δ) [see Eq. (84) below]. This term d * is determined by assuming that the Jacobian D W satisfies the conservation law of the phase-space volume element written as where V (rl) gc and dU/dt are given by Eqs. (46) and (48). We can rewrite Eq. (83) as Noting that the last line of Eq. (84) is of O(δ), we can take the correction term d * as a small quantity of O(δ).
The left-hand side of Eq. (84) represents the derivative of ln D W along the radially local guiding center orbit labeled by the constant parameters (s, W ). Assuming that the guiding center orbit ergodically covers the (θ, ζ, U ) space, D W = (1 + d * )/m is determined by Eq. (84) except for a factor that is an arbitrary function of (s, W ). In order to uniquely specify D W = (1+d * )/m, we impose another constraint, where · · · represents the flux-surface average defined in Eq. (32). Owing to the condition in Eq. (85), d * is given as a small correction and it does not affect the surfaceaveraged velocity integral of the equilibrium distribution function F 0 as shown by where d 3 v is given by Eq. (56) and F 0 is the local Maxwellian defined in Eq. (38) with the equilibrium density n 0 and temperature T 0 given as flux-surface functions. We next define another distribution function f * by and use Eq. (83) to rewrite Eq. (51) in terms of f * as where the differential operator V is defined by The collision term C L (f ) in Eq. (51) is replaced by C L (f * ) in Eq. (88) where the deviation of C L (f * ) from C L (f ) is of O(δ 2 ) and it is neglected. It is shown from Eq. (83) that the differential operator V satisfies the antisymmetry relation, where α and β are arbitrary smooth functions on the phase space.
Replacing f with f * in Eq. (77), we can define modified transport fluxes, Γ ncl * a , q ncl * a , and J * E , the values of which agree with those of Γ ncl a , q ncl a , and J E , respectively, to the lowest order in δ because f * = f [1 + O(δ)]. Then, substituting the solution f * of Eq. (88) into the definitions of the modified transport fluxes, we can derive the neoclassical transport equations relating (Γ ncl * a , q ncl * a , J * E ) to (X b1 , X b2 , X E ). These transport equations take the same forms as those in Eq. (81), and we use (L 11 * ab , L 12 * ab , · · · ) to represent the modified transport coefficients which correspond to (L 11 ab , L 12 ab , · · · ) in Eq. (81), respectively. It is shown in the same way as in Sec. III that no additional sources and/or sinks are required to obtain stationary particle and energy balances from Eq. (88) We now multiply Eq. (88) for particle species a with T a f * a /F a0 and take its velocity-space integral, fluxsurface average, and summation over species. Then, we obtain where the inequality is due to the H-theorem given in Eq. (80). Equation (91) means that the neoclassical transport process is subject to the second law of thermodynamics: the summation of products between the transport fluxes and forces equals the entropy production rate expressed in terms of the linearized collision operator, which is positive definite.
Since the differential operator V and the linearized collision operator C L satisfy the antisymmetry relation in Eq. (90) and the adjointness relations in Eq. (79), respectively, we can use the same procedures as in Ref. 29 to prove that the modified transport coefficients (L 11 * ab , L 12 * ab , · · · ) obey the Onsager symmetry relations written as where β ≡ [E s , ∇s · ∇B, ∇s · (b · ∇b)] represent the variables associated with the perpendicular guiding center velocity (V (rl) gc ) ⊥ as explained after Eq. (81). Note that the change from β to −β corresponds to turning (V (rl) gc ) ⊥ in the opposite direction.
The positive definiteness and the Onsager symmetry shown in Eqs. (91) and (92) for the neoclassical transport defined by the solution f * of Eq. (88) do not hold for the neoclassical fluxes (Γ ncl a , q ncl a , J E ) and the transport coefficients (L 11 ab , L 12 ab , · · · ) in Eq. (81) derived from the solution f of Eq. (51). On the other hand, we also find that Γ ncl * a defined by f * is not written in the same form as in Eq. (70) because the parallel momentum balance equation derived from Eq. (88) cannot be used exactly in the same way as in deriving Eq. (70) from Eq. (65). Thus, the intrinsic ambipolarity condition for axisymmetric and quasi-symmetric systems is slightly broken by the modified neoclassical particle fluxes Γ ncl * a obtained using (L 11 * ab , L 12 * ab , L 1 * aE ) while it is rigorously satisfied by Γ ncl a using (L 11 ab , L 12 ab , L 1 aE ) as shown in Sec. III. By the way, it can be shown in the same way as in Ref. 29 that, for axisymmetric systems with up-down symmetry and helical systems with stellarator symmetry, the neoclassical transport coefficients (L 11 * ab , L 12 * ab , · · · ) satisfy the restricted forms of the Onsager symmetry relations, L ij * ab (β) = L ij * ab (−β) = L ji * ba (β) (i, j = 1, 2), L i * aE (β) = −L i * aE (−β) = L i * Ea (β) (i = 1, 2), L * EE (β) = L * EE (−β). (93)

V. CONCLUSIONS
In this paper, a novel radially local approximation of the drift kinetic equation is presented. The approximated guiding center equations, which are shown in Eq. (46), have no radial drift velocity component but they maintain the E × B drift and the component of the magnetic drift tangential to the flux surface. In addition, they conserve the particle kinetic energy at the expense of the conservation of the magnetic moment. Under this approximation, a new drift kinetic equation is given by Eq. (51) in the conservative form, which has favorable properties for numerical simulation that any additional terms for particle and energy sources are unnecessary for obtaining stationary solutions. Also, it is shown to satisfy the intrinsic ambipolarity condition for neoclassical particle fluxes in axisymmetric and quasi-symmetric toroidal systems. Another radially local drift kinetic equation is presented in Eq. (88), the solution of which equals that of Eq. (46) to the leading order in the expansion with respect to the drift ordering parameter δ defined by the ratio of the gyroradius to the equilibrium scale length. The positive definiteness of the entropy production due to the neoclassical transport fluxes and the Onsager symmetry of the neoclassical transport coefficients are rigorously guaranteed by the solution of Eq. (88) although it does not exactly assure the intrinsic ambipolarity condition for neoclassical particle fluxes in axisymmetric and quasi-symmetric systems. Thus, Eqs. (51) and (88) each have favorable properties which are weakly broken in the other equation. To the lowest order in δ, the neoclassical transport fluxes derived from both solutions of Eqs. (51) and (88) have the same values as each other, and no additional sources and/or sinks are required for those solutions to satisfy stationary particle and energy balances consistently. Therefore, both drift kinetic equations are considered to be practically useful for numerically evaluating the neoclassical transport fluxes with including effects of the E × B and magnetic drift motions tangential to the flux surface in the framework of the radially local approximation. Numerical applications of the present local model are in progress and their results will be reported elsewhere.

ACKNOWLEDGMENTS
This work was supported in part by NIFS/NINS under the Project of Formation of International Network for Scientific Collaborations and in part by the NIFS Collaborative Research Programs (NIFS14KNTT026).