3.5 Flux conservative representation

3.5.1 General formulation

The starting point of the flux conservative representation is the conservation of the total number of particles in the plasma,

    ∫     ∫          3   3
N =    ⋅⋅⋅  f (X, P )d Xd  P
(3.120)

where X and P are respectively coordinates in configuration and momentum spaces. According to the systems which are used in the calculations, X =(ψ, θ,ϕ) and P =(p,ξ,φ),

{   3    ---Rr--
   d X = |∇ψ||ψ^⋅^r|dψd θdϕ
   d3P =p2dpd ξdφ
(3.121)

as shown in Appendix A, one obtains

    ∫     ∫
                            ---Rr------2
N =    ⋅⋅⋅  f (ψ, θ,ϕ,p,ξ,φ)     ||^   ||p dψdθdϕdpd ξdφ
                            |∇ ψ||ψ ⋅^r|
(3.122)

Using the transformation

-1---ξ2-   1--ξ20-
B (ψ, θ) = B0 (ψ )
(3.123)

that results from conservation of the magnetic moment and energy,

     ∫    ∫                     Rr     2 ξ0B (ψ,θ)
N  =   ⋅⋅⋅   f (ψ,θ,ϕ, p,ξ,φ )----||----||p  ξ-B--(ψ-) dψd θdϕdpdξ0dφ
                            |∇ ψ||^ψ ⋅ ^r|     0
(3.124)

because ξdξ∕B(ψ,θ) = ξ00∕B0(ψ). Since at the zero order, f0 is constant along a magnetic field line, f0 = f0(0) is independent of the poloidal angle θ, where f0(0) is the bounce averaged distribution function. Hence

       ∫∫ ∫
             (0)         2
N  =        f0  (ψ,p,ξ0)p dψdpd ξ0 ×
       ∫ θmax                       ∫ 2π    ∫ 2π
             ----Rr|---|-ξ0-B-(ψ,θ)dθ     dϕ     dφ
        θmin  |∇ψ |||ψ^⋅^r|| ξ B0 (ψ )    0      0
          ∫ ∫∫
   =   4π2      f(0)(ψ, p,ξ0) --1---p2dψdpd ξ0 ×
                 0          B0 (ψ )
       ∫ θmax    Rr     ξ0
             -----||---||---B (ψ, θ)dθ                        (3.125)
        θmin  |∇ψ ||ψ^⋅^r| ξ
taking into account, in addition, of the toroidal aximmetry, and cylindrical symmetry of the distribution along the magnetic field line direction. Here, θmin and θmax depends of the particle trajectory in the configuration space, wether they are passing or trapped.

From plasma equilibrium, since

|∇ ψ|
----- = BP
  R
(3.126)

where BP is the poloidal magnetic field,

∫ θmax    Rr     ξ0             ∫ θmax   r   ξ0 B
      ----||----||--B (ψ, θ)dθ =       ||---||------dθ = 2πλ (ψ,ξ0)^q(ψ )
 θmin  |∇ ψ ||ψ^⋅^r| ξ              θmin  |ψ^⋅^r| ξ BP
(3.127)

Here, appears, as expected the normalized bounce time λ(ψ,ξ0) and the factor q^ (ψ) introduced in Sec. 2.2.1, which in conjunction with B0(ψ ) characterizes the local shape of magnetic flux surface. Hence,

        ∫ ∫∫              ^q (ψ )
N  = 8π3     f (00)(ψ,p,ξ0) ------λ(ψ,ξ0)p2dψdpd ξ0
                          B0(ψ )
(3.128)

From this expression, the Jacobian J of the coordinate system (ψ, p,ξ0) may be simply defined as,

               q^(ψ)-         2
J =  JψJpJξ0 = B0 (ψ )λ(ψ, ξ0) p
(3.129)

where

({  Jψ = ^q(ψ )∕B0 (ψ )
   J  = p2
(   p
   Jξ0 = λ (ψ,ξ0)
(3.130)

and the generic conservative form of the kinetic equation may be immediatly deduced

  (0)
∂f0--+ ∇       ⋅S(0) = 0
 ∂t      (ψ,p,ξ0)
(3.131)

where the phase space flux S(0) at B = Bmin is decomposed into a diffusive and a convective part

S (0) = - D (0)∇f (0)+ F(0)f(0)
               0        0
(3.132)

in the mean field theory. Here, D(0) and F(0) are respectively the diffusion tensor and convection vector in phase space. They can be expressed generally as

      (  D (0)  D (0)  D(0) )
 (0)  |    ψ(0ψ)    ψ(0p)   ψ(ξ0) |
D   = (  D pψ   D pp   Dpξ  )
         D (0)  D (0)  D(0)
           ξψ     ξp    ξξ
(3.133)

and

      (   (0) )
 (0)  |  Fp(0) |
F   = (  Fξ   )
         F(0)
          ψ
(3.134)

where each element is function of (ψ, p,ξ0). Here the gradient vector in the reduced (ψ, p,ξ0) space is

     ( ∇   = |∇ψ |∂∕∂ψ     )
     |   ψ                 |
∇ =  ( ∇p  = ∂∕∂√p---2      )
       ∇ ξ0 = ---1-pξ0∂∕∂ξ0
(3.135)

so, following calculations given in Appendix A,

             (    )          (         )        (         )        (         )
∇       ⋅S (0) f (0)   =   1-∂-- JS(0) ⋅ep + 1--∂-- JS(0) ⋅eξ + 1--∂- J S(0) ⋅eψ
  (ψ,p,ξ0)       0         J ∂p              J ∂ξ0              J ∂ψ
                         1 ∂ (   (0)  )   1 1 ∂  ( ∘ ------  (0)  )
                     =   J-∂p- JS   ⋅ ^p - J-p∂-ξ-    1- ξ20JS   ⋅ξ^
                               (              ) 0
                        + 1--∂-  J |∇ ψ| S(0) ⋅ ^ψ
                          J ∂ψ         0              (                   )
                         1  ∂ ( 2 (0))      1    1 ∂   ∘ -----2         (0)
                     =   p2∂p- p Sp   -  λ(ψ,ξ-)-p∂ξ--   1 - ξ0λ(ψ,ξ0)S ξ
                                          (   0     0              )
                        + ---B0-(ψ-)----∂-  -^q(ψ)-λ (ψ, ξ0)|∇ ψ| S (0)       (3.136)
                          λ (ψ,ξ0)^q(ψ )∂ψ   B0 (ψ)            0  ψ
where |∇ ψ|0 is taken on the magnetic flux surface where B is minimum, i.e., B = B0. The first two terms correspond to the usual dynamics in momentum space at a given spatial position ψ, while the third one is associated to spatial transport at fixed p and ξ0. It is interesting to note that spatial transport is not independent of the momentum dynamics through the parameter λ(ψ,ξ0). It corrects the spatial transport from the particle dynamics along the magnetic field line, since most particles tend to spend more time far from B = B0. In the limit of strongly passing particles, λ(ψ,ξ0) 1, and the spatial term becomes independent of ξ0.

It is interesting to cross-check the conservative nature of the transport equation is well ensured by performing the integral

     [                     ]
∫∫ ∫   ∂f(00)            (0)
       -∂t--+ ∇ (ψ,p,ξ0) ⋅S   Jdpd ξ0dψ = 0
(3.137)

or

      ∫ ∫∫
∂N--                  (0)
 ∂t +       ∇(ψ,p,ξ0) ⋅S  Jdpd ξ0dψ =  0
(3.138)

Indeed, in that case, the variation of the total number of particles ∂N∕∂t as a function of time depends only from boundary conditions.

∫∫∫
     ∇ (ψ,p,ξ0) ⋅S(0)J dpdξ0dψ = Ip - Iξ0 + Iψ
(3.139)

where

     ∫∫ ∫       (     )
Ip =      -1 ∂-- p2S(0) Jdpd ξ0dψ
          p2 ∂p     p
(3.140)

     ∫ ∫∫     1    1 ∂  (∘ -----2         (0))
Iξ0 =      λ(ψ,-ξ)-p∂ξ--   1 - ξ0λ (ψ, ξ0)Sξ   Jdpd ξ0d ψ
                 0    0
(3.141)

     ∫ ∫∫     B  (ψ )    ∂ (  ^q(ψ)                 )
Iψ =       ----0----------  ------λ (ψ,ξ0)|∇ψ |S(ψ0)  Jdpd ξ0d ψ
           λ(ψ,ξ0)^q (ψ )∂ψ   B0 (ψ)
(3.142)

For Ip,

   ∫ ∫∫
         1--∂-( 2  (0))
         p2∂p  p S p  J dpdξ0dψ
   ∫ ∫ [      ]pmax
=       p2S(p0)     q^(ψ)-λ(ψ, ξ0) dξ0d ψ
        ∫∫     0   B0 (ψ )
    2        (0)       -^q(ψ)-
=  pmax     Sp (pmax) B0 (ψ )λ (ψ, ξ0)dξ0dψ             (3.143)
and assuming that limp→∞p2Sp(0) = 0, one finds Ip = 0. This condition is generally well fullfiled, except in strong runaway regimes, where above the Dreicer limit characterized by critical momentum pD, electrons gain energy up to very high energies, that are usually well beyond the domain of integration addressed in numerical calculations for the current drive problem. However, in this case, Ip is given by Sp(0) at pmax, where pmax corresponds to the boundary of the momentum domain of integration. Its conservative nature is well ensured, since ∂N∕∂t only depends of this parameter for p.

The integration of ξ0 leads to

        ∫∫ ∫              ( ∘ ------           )
             ---1----1-∂---       2         (0)
Iξ0  =        λ (ψ,ξ0)p ∂ξ0    1- ξ0λ (ψ,ξ0)Sξ   J dpdξ0dψ
        ∫∫  [∘ ------           ]+1
    =          1- ξ2λ (ψ,ξ0)S(0)   -^q(ψ-)pdpd ψ
                   0         ξ   -1B0 (ψ)
    =   0                                                      (3.144)
which indicates that pitch-angle scattering never contributes to variations of N, and therefore the conservative nature of this term in the transport equation is also well demonstrated.

Finally,

       ∫ ∫∫     B0 (ψ)    ∂  ( ^q (ψ )              (0))
Iψ  =        λ(ψ,-ξ)-^q(ψ)∂-ψ  B--(ψ-)λ(ψ,ξ0)|∇ ψ|Sψ   J dpdξ0dψ
       ∫ ∫ [       0            0    ]ψ
             -^q-(ψ-)--             (0)  a 2
    =        Bmin(ψ )λ(ψ,ξ0)|∇ ψ|Sψ   0 p dpdξ0
        q^(ψ )        ∫∫
    =   ----a--|∇ψ |ψa     λ(ψa,ξ0)S(ψ0)(ψa)p2dpdξ0                 (3.145)
        B0 (ψa )
which only depends of edge values at ψa since |∇ψ |0 = R0BP0 = 0 at ψ = 0 when no particle is injected at the plasma center. Here, BP0 is the poloidal magnetic field where B is minimum. If Sψ(0)(ψa) = 0, Iψ = 0, and the total number of particles is conserved in the discharge.

It is important to notice that the magnetic moment is intrinsically conserved in the equations, in particular for the radial transport part, through the pitch-angle dependence of the normalized bounce time λ. Therefore, spatial transport is valid not only for circulating particles satisfying p∕p1, but also for highly trapped electrons, i.e. when p∕p1.

3.5.2 Dynamics in Momentum Space

Momentum space operator It is possible to recover the general bounce averaged transport equation in momentum space by another independent approach. Here, the momentum space dynamics of the kinetic equation is be expressed in conservative form as a flux divergence that may expressed according to the (A.57) introduced in Appendix A

           1  ∂ (     )
∇p ⋅ Sp = -----i JpSip
          Jp ∂p
(3.146)

where Jp is the momentum space Jacobian associated with the momentum space coordinate system (p,ξ,φ ), described in (A.247). Since the spherical system has the natural symmetry of collisions, the momentum space Jacobian is (A.269)

Jp = p2
(3.147)

so that, taking into account that the kinetic equation is gyroaveraged and therefore the coordinate φ disappears, the following expression for the divergence (A.278) is obtained

          -1-∂-( 2  )   1-∂-( ∘ ----2- )
∇p  ⋅Sp = p2∂p  p Sp  - p ∂ξ    1- ξ S ξ
(3.148)

where by definition

Sp = Sp ^p (3.149)
Sξ = Sp ^
ξ (3.150)

In the mean-field theory, the momentum space fluxes may be expressed as the sum of diffusive and convective parts,

Sp = - Dp ∇p ℧ + Fp℧
(3.151)

with

        (            )
D    =     Dpp  Dp ξ                          (3.152)
  p        D ξp  D ξξ
        (  F  )
Fp   =      p                                 (3.153)
           Fξ

The gradient vector p in the reduced coordinates system (p,ξ) in given by (A.277)

     (  ∇p  )
∇p =    ∇
          ξ
(3.154)

with

p =  ∂
∂p- (3.155)
ξ = -∘ ------
--1---ξ2
   p-∂-
∂ξ (3.156)

so that

               ∘ ------
          ∂f     1 - ξ2    ∂f
Sp = - Dpp---+ --------Dp ξ---+ Fpf
          ∂p       p       ∂ξ
(3.157)

and

          ∂f   ∘1----ξ2    ∂f
Sξ = - D ξp--+ --------D ξξ---+ Fξf
          ∂p       p       ∂ξ
(3.158)

Bounce-averaged operator The bounce averaged operator is

            { 1 ∂  (    )}   { 1 ∂ (∘ ------  )}
{∇p ⋅Sp } =   -2--- p2Sp   -   ----   1 - ξ2Sξ
              p ∂p             p∂ ξ
(3.159)

where the bounce averaging operation is defined in (2.62)

         [  ∑  ]  ∫ θmax
{A } = 1-- 1-            dθ|-1--|-r--B-ξ0A
       λ^q  2 σ   T θmin  2π||^ψ ⋅^r||Rp BP  ξ
(3.160)

and ξ is given along the trajectory by

             ∘ -----------(----2)-
ξ(ψ,θ,ξ0) = σ  1-  Ψ(ψ, θ) 1- ξ0
(3.161)

with

          B (ψ,θ)
Ψ (ψ,θ) = --------
          B0 (ψ)
(3.162)

as shown in Sec.2.2.1.

We find from (3.161) that in momentum space

ξdξ = Ψ ξ0dξ0
(3.163)

and we also get

(    2)     (    2)
 1- ξ   = Ψ  1- ξ0
(3.164)

Then, keeping in mind that |ξ0| = σξ0 is independent of σ, we can transform as follows,

{     (  ------  )}
  1-∂- ∘ 1 - ξ2S
  p∂ξ           ξ
  = 1--
λ^q [   ∑  ]
  1-
  2  σT θminθmax dθ-
2π|-1--|
||^ψ ⋅^r||-r-
RpB--
BPξ0
 ξ1-
p∂--
∂ξ(∘ ------  )
   1 - ξ2Sξ
  = 1--
λp--∂--
∂ σξ01-
^q [     ]
 1-∑
 2
    σT θminθmax dθ-
2π--1---
||^   ||
|ψ ⋅^r|-r-
Rp-B-
BP-σ
Ψ(∘ -----2  )
   1-  ξS ξ
  = 1--
λp-∂--
∂ ξ0∘ ------
  1- ξ20σ1-
^q [  ∑  ]
  1-
  2 σT θminθmax dθ-
2π|-1--|
||^ψ ⋅^r||r--
RpB--
BP√σ--
  Ψ(Sξ)
  = 1--
λp-∂--
∂ ξ0∘ ------
  1- ξ2
      0λσ{        }
  √σ-ξ-S ξ
    Ψξ0 (3.165)

using relation ξdξ = Ψξ00 that is deduced from expression (3.161).

Finally, we can rewrite the equation (3.159) in a conservative form as

            1 ∂  ( 2 (0))   1  ∂  (∘  ----2-  (0))
{∇p ⋅Sp} =  p2∂p- p Sp   -  λp∂-ξ-    1- ξ0λS ξ0
                                0
(3.166)

where the following components are defined

  (0)
Sp  =  {Sp}
(3.167)

and

        {        }
 (0)       σ ξ
Sξ0 = σ   √-Ψξ-Sξ
              0
(3.168)

Here, expression (3.166) is completely equivalent to the momentum transport equation deduced from particle conservation. However, this equivalence may be used only because the bounce-averaged operator is local, and does not depends of ψ.

3.5.3 Dynamics in Configuration Space

Configuration space operator The operator that describes the spatial transport is given by relation

  (0)||                   (                        )
∂f0--|| = ----B0(ψ-)----∂-  -^q(ψ)-λ (ψ, ξ0)|∇ ψ |0S (0ψ)
 ∂t  |r  λ (ψ,ξ0)^q(ψ )∂ψ   B0 (ψ)
(3.169)

and since BP = |∇ ψ|∕R, it may be rewritten in the form

     |
∂f(00)||      B0 (ψ)    ∂  (            BP0 (ψ)         (0))
-∂t--|| = λ-(ψ,-ξ-)^q(ψ)-∂ψ- R0 (ψ )^q(ψ) B--(ψ)-λ (ψ, ξ0)Sψ
     r         0                        0
(3.170)

where R0 and BP0 are taken at the poloidal location where the magnetic field is minimum B = B0. Much in the same way,

 (0)        (0)    (0)    (0) (0)
Sψ   =   - D ψψ∇ ψf0 + F ψ f0
                       (0)
     =   - D (ψ0ψ)|∇ ψ |0 ∂f0-+ Fψ(0)f(00)
                      ∂ψ
            (0)             ∂f (0)    (0) (0)
     =   - D ψψR0 (ψ )BP 0(ψ)--0- + Fψ  f0              (3.171)
                             ∂ ψ
where the diffusion cross-terms D(0), Dψp(0),Dξψ(0) and Dψξ(0) between momentum and configuration spaces have been neglected.

Case of circular concentric flux-surfaces In that case, ^
ψ = ^r , and by definition Dψψ(0) = Drr(0), Fψ(0) = Fr(0), since ψ is here just a label. Therefore, using relation |∇ ψ|0∂∕∂ψ = ∂∕∂r,

S(0) =   - D (0)∇ f (0)+ F (0)f(0)
 ψ          ψψ  ψ 0      ψ  0
            (0)      ∂f(00)     (0) (0)
     =   - D ψψ |∇ ψ|0-∂ψ--+ Fψ  f0
                 (0)
            (0)∂f0--    (0) (0)
     =   - D rr ∂r  + Fr  f0
     =   S(0)                                        (3.172)
          r

Furthermore, since

       r  B0
^q(r) = R-B---
        p  P0
(3.173)

one obtains,

  (0)||circ.                 (               )
∂f0--||     =  -RpBP--0--∂-  R0-rλ(r,ξ0)S(r0)
 ∂t  |r        λ (r,ξ0)r ∂ψ   Rp
                 |∇ψ |    ∂  (              )
           =  -------0------  R0rλ (r,ξ0) S(0r)
              λ (r,ξ0)R0r ∂ ψ(               )
              -----1------∂-             (0)
           =  λ (r,ξ0)R0r ∂r  R0r λ(r,ξ0)Sr               (3.174)

Note that R0 may not be here simplified, since it is a function of r, which corresponds to the toroidal configuration. Dynamics in momentum and configuration spaces are also not decoupled, the normalized bounce time λ(r,ξ0) being on both sides of the radial derivative. Only for strongly circulating electrons,

        (0)|circ.
     ∂f-0-||      -1---∂-(     (0))
|ξli0m|→1   ∂t ||    = rR0 ∂r  rR0S r
           r
(3.175)

since lim|ξ0|1λ(r,ξ0) = 1, and the usual cylindrical conservative expression of the radial transport

        (0)||circ.
      ∂f0--|      1-∂-(  (0))
|lξ0im|→1   ∂t  ||   =  r∂r  rSr
           r
(3.176)

is only found in the case ϵ 1, i.e. when R0 Rp.

3.5.4 Bounce-averaged flux calculation

Zero order term: the Fokker-Planck equation

The bounce averaged Fokker-Planck equation is is given in the conservative form by relation (3.166), with the bounce-averaged fluxes (3.167) and (3.168)

Sp(0) = {Sp} (3.177)
Sξ0(0) = σ{   σξ    }
  √----S ξ
    Ψξ0 (3.178)

Because f0 is constant along a magnetic field line, we have f0(p,ξ) = f0(0)(p,ξ0) which is independent of θ and σ. Using the following identities

-{       }
     ∂f0-
 Dpp ∂p = -{Dpp}∂f(00)
 ∂p (3.179)
{    ∘ ------    }
  D  --1---ξ2∂f0-
   pξ    p    ∂ξ = ∘ ------
--1--ξ20-
   p{          }
  √σξ--D
   Ψ ξ0  pξ  (0)
∂f0--
σ∂ξ0 (3.180)
{Fpf0} = {Fp}f0(0) (3.181)
- σ{  σξ      ∂f }
  √----D ξp--0-
   Ψ ξ0    ∂p = -σ{  σξ     }
  √----D ξp
   Ψ ξ0∂f(0)
--0--
 ∂p (3.182)
σ{            ------    }
    σξ     ∘ 1 - ξ2∂f0
  √-----Dξξ------------
    Ψ ξ0       p     ∂ξ =   ------
∘ 1- ξ2
------0-
   pσ{  ξ2    }
  ---2Dξξ
  Ψ ξ0∂f (0)
---0-
σ ∂ξ0 (3.183)
σ{           }
  --σξ--
  √ Ψξ  Fξf0
      0 = σ{         }
  -σ-ξ--
  √ Ψξ F ξ
      0f0(0) (3.184)

we can rewrite

 (0)     (⊬)       (⊬)    (⊬) (⊬)
Sp  = - Dp  ⋅∇ ,ξ⊬℧ ⊬ +  Fp ℧ ⊬
(3.185)

where the bounce averaged flux is decomposed into

      (   (0))
S(p0)=    Sp
         S(ξ00)
(3.186)

with

Sp(0) = -D pp(0)∂f(00)
 ∂p + ∘ ------
--1---ξ20
    pD(0)∂f(00)
 ∂ξ
   0 + Fp(0)f 0(0) (3.187)
Sξ0(0) = -D ξp(0)  (0)
∂f0--
 ∂p + ∘ ------
--1---ξ20
    pDξξ(0)  (0)
∂f0--
 ∂ξ0 + Fξ(0)f 0(0) (3.188)

by defining the diffusion components

Dpp(0) = {Dpp} (3.189)
D(0) = σ{  σ ξ     }
  √----Dp ξ
    Ψξ0 (3.190)
Dξp(0) = σ{          }
  -σ-ξ--
  √ Ψξ D ξp
      0 (3.191)
Dξξ(0) = {  2     }
  ξ--D ξξ
  Ψξ20 (3.192)

and the convection components

Fp(0) = {Fp } (3.193)
Fξ(0) = σ{         }
  √-σξ--Fξ
    Ψ ξ0 (3.194)

where the gradient vector in the reduced (p,ξ )
   0 momentum space is

       ( ∇p   )
∇p,ξ0 =   ∇
           ξ0
(3.195)

with

p = ∂
∂p- (3.196)
ξ0 =  ∘ ------
---1---ξ20-
    p-∂--
∂ ξ0 (3.197)

3.5.5 Up to first order term: the Drift Kinetic equation

In the first-order drift kinetic equation, the momentum space operator

∇p ⋅Sp (f1)
(3.198)

where the fluxes are expressed as (3.151) may be decomposed as

                     (^)
∇p ⋅Sp (f1) = ∇p ⋅Sp  f  + ∇p ⋅Sp (g)
(3.199)

According to (3.166), we can express the bounce-averaged operator as

{         (  )}
 Jp∇p  ⋅Sp  ^f = ∂
---
∂p(      )
  p2^S(p0)-p
--
λ ∂
----
∂ξ0( ∘ ------  (0))
    1 - ξ20λ^Sξ0 (3.200)
{J ∇  ⋅S  (g)}
  p  p   p = ∂--
∂p(      )
  p2S(0)
     p-p-
λ-∂--
∂ξ0( ∘ ------    )
    1 - ξ2λS(0)
         0  ξ0 (3.201)

where we need to evaluate the bounce-averaged fluxes (3.167) and (3.168) for ^
f and g respectively

S^ p(0) = {   ( ) }
 S   f^
   p (3.202)
S^ ξ0(0) = σ{             }
  -σ-ξ--  ( ^)
  √ Ψξ S ξ f
      0 (3.203)

and

Sp(0) = {Sp (g)} (3.204)
Sξ0(0) = σ{            }
  --σξ--
  √ Ψξ S ξ (g)
      0 (3.205)

Because g is constant along a field line, we have g(p,ξ) = g(0)(p,ξ0)which is independent of θ and σ. Therefore, the fluxes for g have exactly the same expression as for f0 in the zero-order equation described in section. This is why the same notation in (3.201) is used, while the fluxes associated with ^f are noted ^S .

Indeed, ^f has an explicit dependence upon θ, which can be isolated as follows:

^f(ψ,θ,p,ξ) = ξ(ψ,θ,ξ0)^f(0)(ψ, p,ξ)
             Ψ (ψ,θ)ξ0         0
(3.206)

with

              pξ0I (ψ )∂f (0)(ψ,p,ξ0)
^f(0)(ψ, p,ξ0) = -----------0---------
              qeB0 (ψ)     ∂ψ
(3.207)

We can note that ^f (0) is antisymmetric in the trapped region, since f0(0) is symmetric and ξ0 is, of course, antisymmetric. As a result, only σf^ (0) can be taken out of the bounce averaging operator. Taking the bounce-average of each term, we find

{        }
       ∂ ^f
 - Dpp ∂p = -σ{         }
   --ξ-
  σΨ ξ0Dpp∂ ^f(0)
 ∂p (3.208)
{    ∘ ------   }
     --1---ξ2∂f^
 Dp ξ   p    ∂ξ = ∘ ------
--1---ξ20
    p{           }
  --ξ2---
  Ψ3 ∕2ξ2 Dpξ
       0∂f^(0)
 ∂ξ0
+ ∘ ------
--1---ξ20
   pσ{            }
    Ψ---1--
  σ Ψ3∕2ξ3Dpξ
         0 ^
f (0) (3.209)
{    }
 Fpf^ = σ{        }
 σ -ξ-Fp
   Ψξ0^f (0) (3.210)
- σ{   σξ     ∂ ^f}
  √-----Dξp---
    Ψ ξ0   ∂p = -{   ξ2     }
  -3∕2-2D ξp
  Ψ   ξ0∂ ^f(0)
-----
 ∂p (3.211)
σ{          ∘ ----2- ^ }
  √σ-ξ-D ξξ--1--ξ--∂f-
    Ψξ0       p    ∂ξ = ∘ -----2
--1---ξ0
    pσ{   3     }
 -σξ-3D ξξ
 Ψ2 ξ0 ^(0)
∂f---
∂ ξ0
+ ∘1----ξ2
-------0
   p{ ξ(Ψ - 1)    }
  ---2-4--D ξξ
    Ψ ξ0^f (0) (3.212)
σ{          }
  -σξ---  ^
  √Ψ-ξ Fξf
      0 = {          }
  --ξ2---
  Ψ3 ∕2ξ2F ξ
       0^
f (0) (3.213)

where the following relation is used

-∂--
∂ξ0-ξ
ξ0 = σ-∂---
∂|ξ0|∘ -----(--------)-
  1- Ψ   1- |ξ |2
--------------0---
       |ξ0|
= σ  2    2
Ψξ0---ξ-
 |ξ|ξ20
= σΨ---1-
|ξ|ξ2
    0 (3.214)

We can therefore rewrite

    (   )
^S(0) ^f(0) =  - ^D (0)⋅∇   ^f(0) + ^F(0)^f(0)
 p              p    p,ξ0       p

where the bounce averaged flux is decomposed into

      (   (0))
S(p0)=    Sp(0)
         Sξ0
(3.215)

with

^
S p(0) = -^
D pp(0)∂ ^f(0)
 ∂p + ∘ ------
--1---ξ20
    p ^
D(0)∂f^(0)
 ∂ξ0 +  ^
Fp(0) ^
f (0) (3.216)
^S ξ0(0) = -^D ξp(0) ^(0)
∂f---
 ∂p + ∘ -----2
--1---ξ0
    pD^ξξ(0)  ^(0)
∂f---
 ∂ξ0 + F^ξ(0)^f (0) (3.217)

by defining the diffusion components

^Dpp(0) = σ{         }
 σ -ξ-Dpp
   Ψξ0 (3.218)
^D(0) = {   ξ2      }
  --3∕2-2Dp ξ
  Ψ   ξ0 (3.219)
^
Dξp(0) = {           }
  --ξ2---
  Ψ3∕2ξ20D ξp (3.220)
^Dξξ(0) = σ{   3    }
 -σξ-3D ξξ
 Ψ2 ξ0 (3.221)

and the convection components

F^p(0) = σ{   ξ    }
 σ ---Fp
   Ψξ0 +   ------
∘ 1 - ξ2
----3--0
  pξ0{ (Ψ - 1)   }
  ---3∕2--Dpξ
   Ψ (3.222)
 ^
Fξ(0) = {          }
  --ξ2---
  Ψ3 ∕2ξ2F ξ
       0 + ∘ ------
--1---ξ20
  pξ3
    0σ{              }
  σξ-(Ψ---1)
     ξ0Ψ2   D ξξ (3.223)

where we use the fact that σξ03 may be taken out of the bounce averaged operator, since ξ03 is an odd function of ξ0. The gradient vector in the reduced (p,ξ0) momentum space is

       (      )
∇p,ξ =   ∇p
   0     ∇ ξ0
(3.224)

with

p = ∂--
∂p (3.225)
ξ0 =    ------
- ∘ 1 - ξ2
--------0-
    p ∂
----
∂ ξ0 (3.226)