Appendix F
Alternative discrete cross-derivatives coefficients

In Sec. 5.4.1, the discretization scheme is chosen so that cross-derivative terms      |
∂2f(00)||
∂p∂ξ0|l+12,i+12,j+12 are all identical and defined at the center of each cell. Consequently, cross-derivatives are determined in a completely symmetric way with respect to neighboring points, which has be proven to be an efficient procedure for most simulation cases [?]. However, the downside of this approach is that internal boundary conditions are no longer satisfied for cross terms, and must be enforced. This is not a satisfactory situation particularly when a large quasilinear diffusion takes place close to these locations of the phase space, in particular at ξ0 = ±1.

Alternative discretization schemes may be considered in order to avoid this problem and keep the treatment of internal boundaries in a uniform way for all derivatives. In this case, D(0) and Dξp(0) are not separated from the derivatives of the distribution function as done in Sec. 5.4.1, so that both ∂∂p(    (0)∂f(0))
  pD pξ -∂0ξ0- and ∂∂ξ0-(∘ ------   (0)∂f(0))
   1 - ξ20λD ξp-∂0p- must directly discretized on the flux grid, according to the general rules. Coefficients are therefore

 ^q          ||(k+1)             ^q      m∑=8
---p2∇p ⋅S(p0)||               = --l+1∕2--    T[m]
B0           l+1∕2,i+1∕2,j+1∕2   B0,l+1∕2 m=1
(F.1)

with

                              (0)||(k+1)
      - p2i+1D (p0p),l+1∕2,i+1,j+1∕2 ∂f∂0p-||            + p2i+1F (p0,l)+1∕2,i+1,j+1∕2f (00,)l+(k1+∕12),i+1,j+1∕2
T[1] =---------------------------l+1∕2,i+1,j+1∕2------------------------------------
                                        Δpi+1∕2
(F.2)

          (0)           ∂f(0)||(k+1)            (0)          (0)(k+1)
      p2iD pp,l+1∕2,i,j+1∕2 -0∂p-||          - p2iFp,l+1∕2,i,j+1∕2f0,l+1∕2,i,j+1∕2
T[2] = ---------------------l+1∕2,i,j+1∕2------------------------------
                                 Δpi+1∕2
(F.3)

        ∘ -----------
          1 - ξ2                            (0)||(k+1)
T [3] = +-------0,j+1∕2-pi+1D (0)            ∂f0--||
           Δpi+1∕2        pξ,l+1∕2,i+1,j+1∕2  ∂ξ0 |l+1∕2,i+1,j+1∕2
(F.4)

          -----------
        ∘ 1 - ξ2                        (0)||(k+1)
T [4] = --------0,j+1∕2pD (0)          ∂f0--|
           Δpi+1∕2    i pξ,l+1∕2,i,j+1∕2  ∂ξ0 ||
                                           l+1∕2,i,j+1∕2
(F.5)

                      ⌊                 (         )            (0)||(k+1)
                      | D(ξ0ξ),l+1∕2,i+1∕2,j+1  1 - ξ20,j+1 λl+1∕2,j+1 ∂f∂ξ00-||
T [5]  =  - --pi+1∕2---|| ------------------------------------------l+1∕2,i+1∕2,j+1-
           λl+1∕2,j+1∕2 |⌈                    .pi+1∕2Δ ξ0,j+1∕2

                 ∘ ---------
            p      1-  ξ2   λl+1∕2,j+1F (0)           f(0)(k+1)       ⌋
         +  -i+1∕2------0,j+1-----------ξ,l+1-∕2,i+1∕2,j+1-0,l+1∕2,i+1∕2,j+1⌉         (F.6)
                                  Δ ξ0,j+1∕2
                    ⌊               (       )          (0)||(k+1)
                    | D (0ξξ),l+1∕2,i+1∕2,j 1 - ξ20,j  λl+1∕2,j ∂f∂0ξ0-||
T [6] =   ---pi+1∕2---|| ------------------------------------l+1∕2,i+1∕2,j-
         λl+1∕2,j+1∕2|⌈                 pi+1∕2Δ ξ0,j+1∕2

                 ∘ -------
           p       1-  ξ2 λl+1∕2,jF (0)         f(0)(k+1)     ⌋
         + -i+1∕2------0,j--------ξ,l+1∕2,i+1∕2,j-0,l+1∕2,i+1∕2,j⌉            (F.7)
                              Δ ξ0,j+1∕2
                 ∘ -----2---                            (0)|(k+1)
 [7]   ---pi+1∕2-----1---ξ0,j+1- l+1∕2,j+1  (0)            ∂f-0-||
T   = λl+1∕2,j+1∕2 Δ ξ0,j+1∕2 λ        Dξp,l+1∕2,i+1∕2,j+1   ∂p ||
                                                           l+1 ∕2,i+1∕2,j+1
(F.8)

                   ∘ -------                          |
  [8]      pi+1∕2     1-  ξ02,j l+1 ∕2,j (0)          ∂f0(0)||(k+1)
T   = - -l+1∕2,j+1∕2-Δ-ξ------λ     D ξp,l+1∕2,i+1∕2,j-∂p--||
        λ             0,j+1∕2                            l+1∕2,i+1∕2,j
(F.9)

There are two methods for calculating derivatives of f0(0)  that appear in the terms T[3],T[4],T[7] and T[8]. As shown in Fig.*** , it is possible to chose values of the distribution function at all flux grid corners,

   (0)||(k+1)           f(0)(k+1)     - f (0)(k+1)
∂f0--||            =  -0,l+1∕2,i+1,j+1---0,l+1∕2,i,j+1-
 ∂p  |                        Δpi+1 ∕2
      l+1∕2,i+1∕2,j+1
(F.10)

   (0)||(k+1)          (0)(k+1)       (0)(k+1)
∂f0--||          =  f0,l+1∕2,i+1,j --f0,l+1∕2,i,j
 ∂p  |                    Δpi+1 ∕2
      l+1∕2,i+1∕2,j
(F.11)

   (0)||(k+1)            (0)(k+1)         (0)(k+1)
∂f0--||            =  f0,l+1∕2,i+1,j+1 --f0,l+1∕2,i+1,j
 ∂ξ0 |                        Δ ξ0,j+1∕2
      l+1∕2,i+1,j+1∕2
(F.12)

   (0)||(k+1)         f(0)(k+1)   - f (0)(k+1)
∂f0--||          =  -0,l+1∕2,i,j+1---0,l+1∕2,i,j-
 ∂ξ0 |                    Δξ0,j+1∕2
      l+1∕2,i,j+1∕2
(F.13)

or to do it in an intermediate way,

   (0)||(k+1 )           (0)(k+1)           (0)(k+1)
∂f-0-||             = f0,l+1∕2,i+3∕2,j+1 --f0,l+1∕2,i-1∕2,j+1
  ∂p |                         Δpi+1 + Δpi
      l+1 ∕2,i+1∕2,j+1
(F.14)

  (0)||(k+1)         (0)(k+1)        (0)(k+1)
∂f0-||           = f0,l+1∕2,i+3∕2,j --f0,l+1∕2,i-1∕2,j
 ∂p |                     Δpi+1 + Δpi
     l+1∕2,i+1∕2,j
(F.15)

   (0)||(k+1 )          f(0)(k+1)       - f (0)(k+1)
∂f-0-||             = -0,l+1∕2,i+1,j+3∕2---0,l+1∕2,i+1,j-1∕2
 ∂ξ0 |                        Δ ξ0,j+1 + Δ ξ0,j
      l+1 ∕2,i+1,j+1∕2
(F.16)

  (0)||(k+1)         (0)(k+1)        (0)(k+1)
∂f0-||           = f0,l+1∕2,i,j+3∕2 --f0,l+1∕2,i,j-1∕2
∂ ξ0 |                    Δ ξ0,j+1 + Δ ξ0,j
     l+1∕2,i,j+1 ∕2
(F.17)

involving flux grids only in the direction that is not considered by the derivative itself. The advantage of this method is that the total number of interpolating coefficients δ is reduced by a factor 2 as compared to the previous method, which represents a significant numerical saving for the calculations. Nevertheless both methods are here detailed.

It is important to notice that this new discretization scheme needs to evaluate now Dξp(0) and D(0) at new grid points, i.e. (l + 1∕2,i+ 1∕2,j + 1), (l + 1 ∕2,i+ 1∕2,j),(l + 1∕2,i,j + 1∕2) and (l + 1∕2,i+ 1,j + 1 ∕2), while in the scheme discussed in Sec. 5.4.1, these coefficients were determined only at the center of the cells (l + 1∕2,i+ 1∕2,j + 1∕2).

Full flux grid interpolation Since the distribution function is defined on the half grid, while fluxes on the full grid, it is necessary to interpolate, because in derivatives F.10-F.13, values of f0(0) are taken on the full grid. In a general way, whatever the detailed value of the weighting factor δ(0) which will be discussed in the Sec. 5.4.3, one may write for the terms proportional to Dpp(0), Fp(0), Dξξ(0) and Fξ(0) as in Sec. 5.4.1

For terms involving D(0) and Dξp(0),two steps interpolation must be carried out,

  (0)(k+1)           (     (0)         )  (0)(k+1)
f0,l+1∕2,i+1,j+1  =    1 - δp,l+1∕2,i+1,j+1  f0,l+1∕2,i+3 ∕2,j+1
                     (0)          (0)(k+1)
                   +δp,l+1∕2,i+1,j+1f0,l+1∕2,i+1∕2,j+1             (F.18)
while
 (0)(k+1)            (     (0)           )  (0)(k+1)
f0,l+1∕2,i+3∕2,j+1  =    1 - δξ,l+1 ∕2,i+3∕2,j+1  f0,l+1∕2,i+3∕2,j+3∕2
                       (0)            (0)(k+1)
                    + δξ,l+1∕2,i+3∕2,j+1f0,l+1∕2,i+3 ∕2,j+1∕2           (F.19)
and
 (0)(k+1)            (     (0)           )  (0)(k+1)
f0,l+1∕2,i+1∕2,j+1  =    1 - δξ,l+1 ∕2,i+1∕2,j+1  f0,l+1∕2,i+1∕2,j+3∕2
                    + δ(0)           f(0)(k+1)                   (F.20)
                       ξ,l+1∕2,i+1∕2,j+1 0,l+1∕2,i+1 ∕2,j+1∕2
which leads to the result
 (0)(k+1)           (    (0)         ) (     (0)           )  (0)(k+1)
f0,l+1∕2,i+1,j+1  =    1- δp,l+1∕2,i+1,j+1   1 - δξ,l+1∕2,i+3∕2,j+1  f0,l+1∕2,i+3∕2,j+3∕2
                     (    (0)         )  (0)            (0)(k+1)
                   +  1- δp,l+1∕2,i+1,j+1  δξ,l+1∕2,i+3∕2,j+1f0,l+1∕2,i+3∕2,j+1 ∕2
                     (0)         (     (0)           ) (0)(k+1)
                   +δp,l+1∕2,i+1,j+1  1- δξ,l+1∕2,i+1∕2,j+1  f0,l+1∕2,i+1∕2,j+3 ∕2
                     (0)          (0)             (0)(k+1)
                   +δp,l+1∕2,i+1,j+1δξ,l+1∕2,i+1∕2,j+1f0,l+1∕2,i+1∕2,j+1∕2        (F.21)

By performing the appropriate index number transformations (i+ 1 → i) and (j + 1 → j) other interpolations are

                 (              ) (                  )
f0(0,)l+(k1+∕12),i,j+1  =    1- δ(p0,)l+1∕2,i,j+1   1-  δξ(0,)l+1∕2,i+1∕2,j+1  f(00,)l+(k1+∕12),i+1∕2,j+3∕2
                   (    (0)       )  (0)            (0)(k+1)
                 +  1- δp,l+1∕2,i,j+1  δξ,l+1∕2,i+1∕2,j+1f0,l+1∕2,i+1∕2,j+1∕2
                   (0)       (     (0)          )  (0)(k+1)
                 +δp,l+1∕2,i,j+1  1- δξ,l+1∕2,i- 1∕2,j+1  f0,l+1∕2,i-1∕2,j+3∕2
                   (0)        (0)             (0)(k+1)
                 +δp,l+1∕2,i,j+1δξ,l+1∕2,i- 1∕2,j+1f0,l+1∕2,i-1∕2,j+1∕2          (F.22)
                 (              ) (                )
f (0)(k+1)     =    1- δ(0)          1 - δ(0)          f (0)(k+1)
 0,l+1∕2,i+1,j        (  p,l+1∕2,i+1,j  )     ξ,l+1∕2,i+3∕2,j   0,l+1∕2,i+3∕2,j+1∕2
                 +  1- δ(0)         δ(0)         f(0)(k+1)
                        p,l+1∕2(,i+1,j   ξ,l+1∕2,i+3∕2),j 0,l+1∕2,i+3∕2,j- 1∕2
                 +δ(0)         1- δ(0)          f(0)(k+1)
                   p,l+1∕2,i+1,j      ξ,l+1∕2,i+1 ∕2,j  0,l+1∕2,i+1∕2,j+1∕2
                 +δ(0)       δ(0)         f (0)(k+1)                    (F.23)
                   p,l+1∕2,i+1,j ξ,l+1∕2,i+1∕2,j 0,l+1∕2,i+1∕2,j-1∕2
               (            ) (                )
f (0)(k+1)   =    1- δ(0)        1-  δ(0)          f (0)(k+1)
 0,l+1∕2,i,j        (  p,l+1∕2,i,j  )    ξ,l+1∕2,i+1∕2,j  0,l+1∕2,i+1∕2,j+1∕2
               +  1- δ(0)       δ(0)         f(0)(k+1)
                      p,l+1(∕2,i,j   ξ,l+1∕2,i+1∕2),j 0,l+1∕2,i+1 ∕2,j- 1∕2
               +δ(0)      1 - δ(0)          f(0)(k+1)
                 p,l+1∕2,i,j      ξ,l+1∕2,i- 1∕2,j  0,l+1∕2,i- 1∕2,j+1∕2
               +δ(0)     δ(0)         f(0)(k+1)                      (F.24)
                 p,l+1∕2,i,j ξ,l+1∕2,i-1∕2,j 0,l+1∕2,i-1∕2,j-1∕2

Gathering all terms, corresponding matrix coefficients for the zero order Fokker-Planck equation may be expressed as

             |
-^q-p2∇  ⋅S (0)||
B0     p   p |l+12,i+12,j+12(k+1)
=  ^ql+1∕2
--------
B0,l+1∕2 i=i-1i=i+1 j=j-1j=j+1M p,l+12,i+12,j+12(0)f 0,l+12,i+12,j+12(0)(k+1) (F.25)

where

---                  ---                  ---
M (p0,)l+1∕2,i′+1 ∕2,j′+1∕2 = M (0p,)l+[11]∕2,i′+1∕2,j′+1∕2 + M (0p,)l+[21]∕2,i′+1∕2,j′+1∕2

Terms proportional to Dξξ(0), Fξ(0) and Dpp(0), Fp(0) are all gathered in matrix Mp(0)[1] and may be easily obtained from Mp(0) given in Sec. 5.4.1 by taking D(0) = Dξp(0) = 0. The other terms proportional to D(0) and Dξp(0) are

                           ∘1---ξ2------
---(0)[2]                  --------0,j+1∕2--     (0)
M  p,l+1 ∕2,i+3∕2,j+3∕2  =  + Δpi+1 ∕2Δ ξ0,j+1∕2pi+1D pξ,l+1∕2,i+1,j+1∕2
                         (     (0)         )(     (0)           )
                       ×  1 - δp,l+1∕2,i+1,j+1   1- δξ,l+1∕2,i+3∕2,j+1
                                      ∘ ---------
                            pi+1∕2       1 - ξ20,j+1    l+1 ∕2,j+1  (0)
                       + λl+1∕2,j+1∕2Δ-ξ-----Δp------λ       D ξp,l+1∕2,i+1 ∕2,j+1
                         (             0,j+1∕2)(  i+1∕2            )
                       ×  1 - δ(0)            1- δ(0)                    (F.26)
                               p,l+1∕2,i+1,j+1       ξ,l+1∕2,i+3∕2,j+1
                             -----------
                           ∘ 1- ξ2
M--(0)[2]             =  + --------0,j+1∕2--p   D(0)
   p,l+1 ∕2,i+1∕2,j+3∕2       Δpi+1 ∕2Δ ξ0,j+1∕2 i+1  pξ,l+1∕2,i+1,j+1∕2
                          (0)          (     (0)           )
                       × δp,l+1∕2,i+1,j+1 1 - δξ,l+1∕2,i+1∕2,j+1
                           ∘ ----2------
                             1- ξ0,j+1∕2     (0)
                       - Δp-----Δ-ξ------piDpξ,l+1∕2,i,j+1∕2
                         (  i+1 ∕2   0,j+1∕2)(                   )
                       ×  1 - δ(0)          1- δ(0)
                               p,l+1∕2,i,j∘+1-------ξ,l+1∕2,i+1∕2,j+1
                            p           1 - ξ20,j+1
                       + ----i+1∕2-------------------λl+1 ∕2,j+1D (ξ0p),l+1∕2,i+1 ∕2,j+1
                         λl+1∕2,j+1∕2Δ ξ(0,j+1∕2Δpi+1 ∕2     )
                       × δ(0)           1 - δ(0)
                          p,l+1∕2,i+1,j+1∘ ----ξ,l+1∕2,i+1∕2,j+1
                                        1 - ξ2
                       - ---pi+1∕2------------0,j+1---λl+1 ∕2,j+1D (0)
                         λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2          ξp,l+1∕2,i+1 ∕2,j+1
                         (     (0)       )(     (0)           )
                       ×  1 - δp,l+1∕2,i,j+1   1- δξ,l+1∕2,i+1∕2,j+1            (F.27)
                           ∘ -----------
---                          1- ξ2
M--(0)[2]             =  - --------0,j+1∕2--p D(0)
   p,l+1 ∕2,i-1∕2,j+3∕2       Δpi+1 ∕2Δ ξ0,j+1∕2 i pξ,l+1∕2,i,j+1∕2
                          (0)        (     (0)           )
                       × δp,l+1∕2,i,j+1  1 - δξ,l+1 ∕2,i-1∕2,j+1
                                      ∘ -----2---
                         ---pi+1∕2-------1---ξ0,j+1--- l+1 ∕2,j+1  (0)
                       - λl+1∕2,j+1∕2Δ ξ     Δp      λ       D ξp,l+1∕2,i+1 ∕2,j+1
                                    (  0,j+1∕2   i+1∕2   )
                       × δ(p0),l+1∕2,i,j+1  1 - δ(0ξ),l+1 ∕2,i-1∕2,j+1                  (F.28)
                           ∘1---ξ2------
---(0)[2]                  --------0,j+1∕2--     (0)
M  p,l+1 ∕2,i+3∕2,j+1∕2  =  + Δpi+1 ∕2Δ ξ0,j+1∕2pi+1D pξ,l+1∕2,i+1,j+1∕2
                         (     (0)         ) (0)
                       ×  1 - δp,l+1∕2,i+1,j+1 δξ,l+1∕2,i+3∕2,j+1
                           ∘ -----------
                             1- ξ20,j+1∕2       (0)
                       - Δp-----Δ-ξ------pi+1D pξ,l+1∕2,i+1,j+1∕2
                         (  i+1 ∕2   0,j+1∕2)(                 )
                       ×  1 - δ(0)          1- δ(0)
                               p,l+1∕2,i+1∘,j-------ξ,l+1∕2,i+3∕2,j
                            p           1 - ξ20,j+1
                       + ----i+1∕2-------------------λl+1 ∕2,j+1D (ξ0p),l+1∕2,i+1 ∕2,j+1
                         λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2
                         (     (0)         ) (0)
                       ×  1 - δp,l+1∕2,i+1,j+1-δξ,l+1∕2,i+3∕2,j+1
                                       ∘      2
                       - ---pi+1∕2--------1---ξ0,j----λl+1 ∕2,jD (0)
                         λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2        ξp,l+1∕2,i+1∕2,j
                         (     (0)       )(     (0)         )
                       ×  1 - δp,l+1∕2,i+1,j   1- δξ,l+1∕2,i+3∕2,j             (F.29)
                     ∘1----ξ2-----
--(0)[2]              --------0,j+1∕2--     (0)
M p,i+1∕2,j+1∕2 =   + Δpi+1∕2Δ ξ0,j+1∕2 pi+1D pξ,l+1∕2,i+1,j+1∕2
                     (0)          (0)
                  × δp,l+1∕2,i+1,j+1δξ,l+1∕2,i+1∕2,j+1
                     ∘ -----2-----
                    ---1---ξ0,j+1∕2--     (0)
                  - Δpi+1∕2Δ ξ0,j+1∕2 pi+1D pξ,l+1∕2,i+1,j+1∕2
                               (                )
                  × δp(0,)l+1∕2,i+1,j  1- δ(ξ0,)l+1∕2,i+1∕2,j
                     ∘ -----------
                       1 - ξ20,j+1∕2
                  - ----------------piD (p0)ξ,l+1∕2,i,j+1∕2
                    Δ(pi+1∕2Δ ξ0,j+1∕2)
                  ×  1 - δ(0)         δ(0)
                     ∘ ---p,l+1-∕2,i,j+1   ξ,l+1∕2,i+1∕2,j+1
                       1 - ξ2
                  + --------0,j+1∕2--piD (0)
                    Δpi+1∕2Δ ξ0,j+1∕2    pξ,l+1∕2,i,j+1∕2
                    (     (0)     ) (     (0)         )
                  ×  1 - δp,l+1 ∕2,i,j  1 - δξ,l+1 ∕2,i+1∕2,j
                                 ∘ -----2---
                    --pi+1∕2-------1---ξ0,j+1--- l+1∕2,j+1  (0)
                  + λl+1∕2,j+1∕2 Δξ0,j+1∕2Δpi+1∕2λ        D ξp,l+1 ∕2,i+1∕2,j+1
                     (0)          (0)
                  × δp,l+1∕2,i+1,j+1δξ,l+1∕2,i+1∕2,j+1
                                 ∘ ---------
                      pi+1∕2       1 - ξ20,j+1    l+1∕2,j+1  (0)
                  - λl+1∕2,j+1∕2-Δξ------Δp-----λ        D ξp,l+1 ∕2,i+1∕2,j+1
                    (            0,j+)1∕2   i+1∕2
                  ×  1 - δ(0p,)l+1 ∕2,i,j+1  δ(ξ0,l)+1∕2,i+1∕2,j+1
                                  ∘ -------
                      p             1 - ξ20,j
                  - ---i+1∕2------------------λl+1∕2,jD (ξ0)p,l+1∕2,i+1∕2,j
                    λl+1∕2,j+1∕2 Δ(ξ0,j+1∕2Δpi+1∕2  )
                  × δ(0)         1- δ(0)
                    p,l+1∕2,i+1,j    ∘ ξ,l+1∕2,i+1∕2,j
                                    1 - ξ2
                  + --pi+1∕2-------------0,j---λl+1∕2,jD (0)
                    λl+1∕2,j+1∕2 Δξ0,j+1∕2Δpi+1∕2         ξp,l+1∕2,i+1∕2,j
                    (     (0)     ) (     (0)         )
                  ×  1 - δp,l+1 ∕2,i,j  1 - δξ,l+1 ∕2,i+1∕2,j                  (F.30)
                           ∘ ----2------
---(0)[2]                      1- ξ0,j+1∕2     (0)
M  p,l+1 ∕2,i-1∕2,j+1∕2  =  - Δp-----Δ-ξ------piDpξ,l+1∕2,i,j+1∕2
                            i+1 ∕2   0,j+1∕2
                       × δ(p0),l+1∕2,i,j+1 δ(0ξ,)l+1∕2,i-1∕2,j+1
                           ∘ -----------
                             1- ξ20,j+1∕2     (0)
                       + ----------------piDpξ,l+1∕2,i,j+1∕2
                         Δpi+1 ∕2Δ ξ(0,j+1∕2          )
                       × δ(0)       1 - δ(0)
                          p,l+1∕2,i,j    ∘ ξ,l+1∕2,i-1∕2,j
                                        1 - ξ2
                       - ---pi+1∕2------------0,j+1---λl+1 ∕2,j+1D (0)
                         λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2          ξp,l+1∕2,i+1 ∕2,j+1
                          (0)         (0)
                       × δp,l+1∕2,i,j+1 δξ,l+1∘∕2,i-1∕2,j+1
                                         1 - ξ2
                       + ---pi+1∕2-------------0,j----λl+1 ∕2,jD (0)
                         λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2        ξp,l+1∕2,i+1∕2,j
                          (0)      (     (0)         )
                       × δp,l+1∕2,i,j 1 - δξ,l+1∕2,i-1∕2,j                      (F.31)
                          ∘ -----------
--(0)[2]                      1- ξ20,j+1∕2
M p,l+1∕2,i+3∕2,j-1∕2  =  - ----------------pi+1Dp(0ξ),l+1∕2,i+1,j+1∕2
                        Δ(pi+1 ∕2Δ ξ0,j+1∕2)
                      ×   1- δ(0)        δ(0)
                              p,l+1∕2,i+1∘,j -ξ,l+1∕2,i+3∕2,j
                                         1- ξ2
                      - ---pi+1∕2-------------0,j----λl+1∕2,jD (0)
                        λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2        ξp,l+1∕2,i+1 ∕2,j
                        (     (0)       ) (0)
                      ×   1- δp,l+1∕2,i+1,j δξ,l+1∕2,i+3∕2,j                (F.32)
                          ∘ -----------
---                         1- ξ2
M-(0)[2]             =  - --------0,j+1∕2--p   D (0)
  p,l+1∕2,i+1∕2,j-1∕2       Δpi+1 ∕2Δ ξ0,j+1∕2 i+1 pξ,l+1∕2,i+1,j+1∕2
                         (0)         (0)
                      × δp,l+1∕2,i+1,jδξ,l+1∕2,i+1∕2,j
                          ∘ ----2------
                        ----1--ξ0,j+1∕2--   (0)
                      + Δpi+1 ∕2Δ ξ0,j+1∕2piDpξ,l+1∕2,i,j+1∕2
                        (     (0)     ) (0)
                      ×   1- δp,l+1∕2,i,j δξ,l+1∕2,i+1∕2,j
                                      ∘  -------
                           pi+1∕2         1- ξ20,j            (0)
                      - -l+1∕2,j+1∕2----------------λl+1∕2,jD ξp,l+1∕2,i+1 ∕2,j
                        λ          Δ ξ0,j+1∕2Δpi+1 ∕2
                      × δ(0)        δ(0)
                         p,l+1∕2,i+1,j ξ,l+1∘∕2,i+1∕2,j-
                           p             1- ξ20,j
                      + -l+1i+∕12∕,j2+1∕2----------------λl+1∕2,jD (ξ0p),l+1∕2,i+1 ∕2,j
                        λ(          Δ ξ0,)j+1∕2Δpi+1 ∕2
                      ×   1- δ(0)      δ(0)                           (F.33)
                              p,l+1∕2,i,j  ξ,l+1∕2,i+1∕2,j
                          ∘1---ξ2------
--(0)[2]                  --------0,j+1∕2--   (0)
M p,l+1∕2,i-1∕2,j-1∕2  =  + Δpi+1 ∕2Δ ξ0,j+1∕2piDpξ,l+1∕2,i,j+1∕2
                         (0)       (0)
                      × δp,l+1∕2,i,jδξ,l+1∕2,i-1∕2,j
                                      ∘  ----2--
                        ---pi+1∕2---------1--ξ0,j---- l+1∕2,j (0)
                      + λl+1∕2,j+1∕2Δ ξ0,j+1∕2Δpi+1 ∕2λ     D ξp,l+1∕2,i+1 ∕2,j
                         (0)       (0)
                      × δp,l+1∕2,i,jδξ,l+1∕2,i-1∕2,j                         (F.34)

Considerable simplifications occur for uniform momentum and pitch-angle grids. In that case, δξ(0) = δp(0) = 12, while Δp and Δξ0 are constant,

                         ∘ -----------
---                        1- ξ2
M-(0)[2u]            =   + ------0,j+1∕2p   D (0)
  p,l+1∕2,i+3∕2,j+3 ∕2          4Δp Δ ξ0    i+1  pξ,l+1∕2,i+1,j+1∕2
                                   ∘ -----2---
                         --pi+1∕2-----1---ξ0,j+1  l+1∕2,j+1  (0)
                       + λl+1∕2,j+1∕2  4Δ ξ0Δp  λ        D ξp,l+1∕2,i+1∕2,j+1(F.35)
                         ∘1---ξ2------
--(0)[2u]                  ------0,j+1∕2-     (0)
M p,l+1∕2,i+1∕2,j+3∕2  =   +   4Δp Δ ξ0  pi+1D pξ,l+1∕2,i+1,j+1∕2
                         ∘ ----2------
                         --1--ξ0,j+1∕2-   (0)
                       -   4Δp Δ ξ0  piD pξ,l+1∕2,i,j+1∕2          (F.36)
                         ∘ -----------
--(0)[2u]                    1- ξ2
M p,l+1∕2,i-1∕2,j+3 ∕2  =   - ------0,j+1∕2piD (0)
                           4Δp Δ ξ0 ∘ ----pξ,l+1∕2,i,j+1∕2
                                     1 - ξ2
                       - --pi+1∕2----------0,j+1λl+1∕2,j+1D (0)           (F.37)
                         λl+1∕2,j+1∕2  4Δ ξ0Δp             ξp,l+1∕2,i+1∕2,j+1
                                   ∘ ---------
--(0)[2u]                    pi+1∕2     1 - ξ20,j+1  l+1∕2,j+1  (0)
M p,l+1∕2,i+3∕2,j+1 ∕2  =   + -l+1∕2,j+1∕2--4Δ-ξ-Δp--λ        D ξp,l+1∕2,i+1∕2,j+1
                         λ         ∘ ----0--
                           p         1 - ξ20,j
                       - -l+1i∕+21,j∕2+1∕2---------λl+1∕2,jD (0ξp),l+1∕2,i+1∕2,j     (F.38)
                         λ          4Δ ξ0Δp
--(0)[2u]
M p,i+1∕2,j+1∕2 = 0

                                   ∘ -----2---
--(0)[2u]                  --pi+1∕2-----1---ξ0,j+1  l+1∕2,j+1  (0)
M p,l+1∕2,i-1∕2,j+1 ∕2  =   - λl+1∕2,j+1∕2  4Δ ξ0Δp  λ        D ξp,l+1∕2,i+1∕2,j+1
                                   ∘ -----2-
                         --pi+1∕2-----1---ξ0,j l+1∕2,j  (0)
                       + λl+1∕2,j+1∕2 4Δ ξ0Δp λ      D ξp,l+1∕2,i+1∕2,j     (F.39)
                         ∘ -----------
--(0)[2u]                    1- ξ2
M p,l+1∕2,i+3∕2,j-1∕2  =   - ------0,j+1∕2pi+1D (0)
                           4Δp Δ ξ0∘  -----pξ,l+1∕2,i+1,j+1∕2
                                      1- ξ2
                       - --pi+1∕2---------0,jλl+1∕2,jD (0)            (F.40)
                         λl+1∕2,j+1∕2 4Δ ξ0Δp          ξp,l+1∕2,i+1∕2,j
                         ∘ -----------
---                        1- ξ2
M-(0)[2u]            =   - ------0,j+1∕2p   D (0)
  p,l+1∕2,i+1∕2,j-1∕2          4Δp Δ ξ0   i+1  pξ,l+1∕2,i+1,j+1∕2
                         ∘ ----2------
                         --1--ξ0,j+1∕2-   (0)
                       +   4Δp Δ ξ0  piD pξ,l+1∕2,i,j+1∕2          (F.41)
                         ∘ -----------
--(0)[2u]                    1- ξ20,j+1∕2
M p,l+1∕2,i-1∕2,j-1∕2  =   + ------------piD (0pξ),l+1∕2,i,j+1∕2
                           4Δp Δ ξ0∘  -------
                                      1- ξ2
                       + --pi+1∕2---------0,jλl+1∕2,jD (0)            (F.42)
                         λl+1∕2,j+1∕2 4Δ ξ0Δp          ξp,l+1∕2,i+1∕2,j

Half flux grid interpolation In this case, a single interpolation step is necessary,

 (0)(k+1)            (     (0)           )  (0)(k+1)
f0,l+1∕2,i+3∕2,j+1  =    1 - δξ,l+1 ∕2,i+3∕2,j+1  f0,l+1∕2,i+3∕2,j+3∕2
                    + δ(0)           f(0)(k+1)                   (F.43)
                       ξ,l+1∕2,i+3∕2,j+1 0,l+1∕2,i+3 ∕2,j+1∕2
 (0)(k+1)            (     (0)           )  (0)(k+1)
f0,l+1∕2,i-1∕2,j+1  =    1 - δξ,l+1 ∕2,i-1∕2,j+1  f0,l+1∕2,i- 1∕2,j+3∕2
                       (0)            (0)(k+1)
                    + δξ,l+1∕2,i-1∕2,j+1f0,l+1∕2,i- 1∕2,j+1∕2           (F.44)
 (0)(k+1)          (     (0)         )  (0)(k+1)
f0,l+1∕2,i+3∕2,j =    1 - δξ,l+1 ∕2,i+3∕2,j  f0,l+1∕2,i+3∕2,j+1∕2
                     (0)          (0)(k+1)
                  + δξ,l+1∕2,i+3∕2,jf0,l+1∕2,i+3∕2,j- 1∕2             (F.45)
 (0)(k+1)          (     (0)         )  (0)(k+1)
f0,l+1∕2,i-1∕2,j =    1 - δξ,l+1 ∕2,i-1∕2,j  f0,l+1∕2,i-1∕2,j+1∕2
                  + δ(0)         f(0)(k+1)                     (F.46)
                     ξ,l+1∕2,i-1∕2,j 0,l+1∕2,i- 1∕2,j- 1∕2
and
 (0)(k+1)            (     (0)           )  (0)(k+1)
f0,l+1∕2,i+1,j+3∕2 =    1 - δp,l+1 ∕2,i+1,j+3 ∕2  f0,l+1∕2,i+3∕2,j+3∕2
                       (0)            (0)(k+1)
                    + δp,l+1∕2,i+1,j+3∕2f0,l+1∕2,i+1 ∕2,j+3∕2           (F.47)
 (0)(k+1)            (     (0)           )  (0)(k+1)
f0,l+1∕2,i+1,j-1∕2 =    1 - δp,l+1 ∕2,i+1,j- 1∕2  f0,l+1∕2,i+3∕2,j-1∕2
                       (0)            (0)(k+1)
                    + δp,l+1∕2,i+1,j-1∕2f0,l+1∕2,i+1 ∕2,j- 1∕2           (F.48)
 (0)(k+1)          (     (0)         )  (0)(k+1)
f0,l+1∕2,i,j+3∕2 =    1 - δp,l+1 ∕2,i,j+3∕2  f0,l+1∕2,i+1∕2,j+3∕2
                  + δ(0)         f(0)(k+1)                     (F.49)
                     p,l+1∕2,i,j+3∕2 0,l+1∕2,i- 1∕2,j+3∕2
 (0)(k+1)          (     (0)         )  (0)(k+1)
f0,l+1∕2,i,j-1∕2 =    1 - δp,l+1 ∕2,i,j- 1∕2  f0,l+1∕2,i+1∕2,j-1∕2
                     (0)          (0)(k+1)
                  + δp,l+1∕2,i,j-1∕2f0,l+1∕2,i- 1∕2,j- 1∕2             (F.50)

Matrix coefficients related to cross-derivatives are then

                              ∘ -----------
--(0)[u]                         1-  ξ02,j+1 ∕2
M p,l+1∕2,i+3∕2,j+3∕2  =   +------------------------pi+1D(p0ξ),l+1∕2,i+1,j+1∕2
                        Δp(i+1 ∕2(Δξ0,j+1 + Δ ξ)0,j)
                       ×  1- δ(0)
                              p,l+1∕2,i+1,j+3∘∕2--------
                                           1-  ξ2
                       +---pi+1∕2---------------0,j+1------λl+1∕2,j+1D (0)
                        λl+1∕2,j+1∕2Δ ξ0,j+1∕2(Δpi+1 + Δpi)           ξp,l+1 ∕2,i+1∕2,j+1
                         (    (0)           )
                       ×  1- δξ,l+1∕2,i+3∕2,j+1                                 (F.51)
                             ∘ -----------
---                            1- ξ2
M--(0)[u]             =  ------------0,j+1∕2------p   D(0)
   p,l+1∕2,i+1∕2,j+3∕2     Δpi+1 ∕2(Δξ0,j+1 + Δ ξ0,j) i+1  pξ,l+1∕2,i+1,j+1∕2
                          (0)
                       × δp,l+1 ∕2,i+1,j+3∕2
                              ∘  ----2------
                         --------1--ξ0,j+1∕2------   (0)
                       - Δpi+1∕2(Δ ξ0,j+1 + Δξ0,j)piD pξ,l+1∕2,i,j+1∕2
                         (     (0)         )
                       ×  1 - δp,l+1∕2,i,j+3∕2                          (F.52)
                              ∘ -----------
--(0)[u]                         1-  ξ2
M p,l+1∕2,i-1∕2,j+3∕2  =   -------------0,j+1-∕2------piD(0)
                        Δpi+1 ∕2(Δξ0,j+1 + Δ ξ0,j)   pξ,l+1∕2,i,j+1∕2
                         (0)
                       ×δp,l+1∕2,i,j+3∕2    ∘ ---------
                                           1-  ξ2
                       ----pi+1∕2---------------0,j+1------λl+1∕2,j+1D (0)
                        λl+1∕2,j+1∕2Δ ξ0,j+1∕2(Δpi+1 + Δpi)           ξp,l+1 ∕2,i+1∕2,j+1
                         (    (0)           )
                       ×  1- δξ,l+1∕2,i-1∕2,j+1                                 (F.53)
                                          ∘ -----2-
--(0)[u]                 ---pi+1∕2------------1---ξ0,j------  l+1∕2,j  (0)
M p,l+1∕2,i+3∕2,j+1∕2  =   -λl+1∕2,j+1∕2Δ ξ0,j+1∕2(Δpi+1 + Δpi)λ      D ξp,l+1∕2,i+1∕2,j
                         (                )
                       ×  1- δ(ξ0,)l+1∕2,i+3∕2,j
                                         ∘ ---------
                           pi+1∕2           1-  ξ02,j+1                 (0)
                       +--l+1∕2,j+1∕2----------------------λl+1∕2,j+1D ξp,l+1 ∕2,i+1∕2,j+1
                        λ          Δ ξ0,j+1∕2(Δpi+1 + Δpi)
                       ×δ(0)                                                 (F.54)
                         ξ,l+1∕2,i+3∕2,j+1
--(0)[u]
M p,i+1∕2,j+1∕2 = 0
(F.55)

                                        ∘ -----2-
--(0)[u]                --pi+1∕2-----------1---ξ0,j------- l+1∕2,j (0)
M p,l+1∕2,i-1∕2,j+1∕2  =   λl+1∕2,j+1 ∕2 Δξ0,j+1∕2(Δpi+1 + Δpi )λ      Dξp,l+1∕2,i+1 ∕2,j
                         (    (0)         )
                       ×  1- δξ,l+1∕2,i-1∕2,j
                                         ∘ ---------
                           pi+1∕2           1-  ξ02,j+1        l+1∕2,j+1  (0)
                       -λl+1∕2,j+1∕2Δ-ξ------(Δp----+-Δp-)λ        D ξp,l+1 ∕2,i+1∕2,j+1
                                      0,j+1∕2    i+1      i
                       ×δ(ξ0,l)+1∕2,i-1∕2,j+1                                       (F.56)
                                          ∘  -------
---                                          1- ξ2
M--(0)[u]             =  - ---pi+1∕2----------------0,j-------λl+1∕2,jD(0)
   p,l+1∕2,i+3∕2,j-1∕2       λl+1∕2,j+1∕2 Δ ξ0,j+1∕2 (Δpi+1 + Δpi )        ξp,l+1∕2,i+1∕2,j
                          (0)
                       × δξ,l+1 ∕2,i+3∕2,j-------
                              ∘  1- ξ2
                       - ------------0,j+1∕2-----p   D (0)
                         Δpi+1∕2(Δ ξ0,j+1 + Δξ0,j) i+1  pξ,l+1 ∕2,i+1,j+1∕2
                         (     (0)           )
                       ×  1 - δp,l+1∕2,i+1,j-1∕2                              (F.57)
                              ∘1---ξ2------
--(0)[u]                  ------------0,j+1∕2------     (0)
M p,l+1∕2,i+1∕2,j-1∕2  =  - Δpi+1 ∕2 (Δ ξ0,j+1 + Δ ξ0,j)pi+1D pξ,l+1∕2,i+1,j+1∕2
                         (0)
                      × δp,l+1∕2,i+1,j-1∕2
                              ∘ ----2------
                        --------1--ξ0,j+1∕2------   (0)
                      + Δpi+1 ∕2 (Δ ξ0,j+1 + Δ ξ0,j)piD pξ,l+1∕2,i,j+1 ∕2
                        (                 )
                      ×   1- δ(p0,l)+1∕2,i,j-1∕2                          (F.58)
                                         ∘ -------
---                                        1- ξ2
M-(0)[u]            =   ---pi+1∕2----------------0,j-------λl+1∕2,jD (0)
  p,l+1∕2,i-1∕2,j-1∕2     λl+1∕2,j+1∕2Δ ξ0,j+1∕2(Δpi+1 + Δpi)         ξp,l+1∕2,i+1∕2,j
                          (0)
                       × δξ,l+1∕2,i-1∕2,j------
                              ∘ 1 - ξ2
                       + ------------0,j+1∕2-----p D (0)
                         Δpi+1∕2(Δ ξ0,j+1 + Δ ξ0,j) i pξ,l+1∕2,i,j+1∕2
                          (0)
                       × δp,l+1∕2,i,j-1∕2                                    (F.59)
and in the special case of uniform momentum and pitch-angle grids

                         ∘1---ξ2------
--(0)[2u]                  ------0,j+1∕2-     (0)
M p,l+1∕2,i+3∕2,j+3 ∕2  =   +   4Δp Δ ξ0   pi+1D pξ,l+1∕2,i+1,j+1∕2
                                   ∘ -----2---
                         --pi+1∕2-----1---ξ0,j+1  l+1∕2,j+1  (0)
                       + λl+1∕2,j+1∕2  4Δ ξ0Δp  λ        D ξp,l+1∕2,i+1∕2,j+1(F.60)
                       ∘ -----------
--(0)[2u]                  1- ξ2
M p,l+1∕2,i+1 ∕2,j+3∕2  =   ------0,j+1∕2pi+1D (0)
                        ∘4Δp-Δ-ξ0-----   pξ,l+1 ∕2,i+1,j+1∕2
                           1- ξ2
                       --------0,j+1∕2piD (0)                   (F.61)
                           4Δp Δξ0      pξ,l+1∕2,i,j+1∕2
                         ∘ -----------
--(0)[2u]                    1- ξ20,j+1∕2    (0)
M p,l+1∕2,i-1∕2,j+3 ∕2  =   - --4Δp-Δ-ξ---piD pξ,l+1∕2,i,j+1∕2
                                 0 ∘ ---------
                           p         1 - ξ20,j+1
                       - -l+1i∕+21,j∕2+1∕2-----------λl+1∕2,j+1D (ξ0)p,l+1∕2,i+1∕2,j+1(F.62)
                         λ           4Δ ξ0Δp
                                   ∘ -------
---                                  1 - ξ2
M-(0)[2u]            =   - --pi+1∕2----------0,jλl+1∕2,jD (0)
  p,l+1∕2,i+3∕2,j+1 ∕2        λl+1∕2,j+1∕2 4Δ ξ0Δp          ξp,l+1∕2,i+1∕2,j
                                   ∘ -----2---
                         --pi+1∕2-----1---ξ0,j+1  l+1∕2,j+1  (0)
                       + λl+1∕2,j+1∕2  4Δ ξ0Δp  λ        D ξp,l+1∕2,i+1∕2,j+1(F.63)
M-(0)[2u]      = 0
  p,i+1∕2,j+1∕2
(F.64)

                                  ∘ -------
---                                 1- ξ2
M (p0,l)[+2u1]∕2,i-1∕2,j+1 ∕2  =   --pi+1∕2---------0,jλl+1∕2,jD (0)
                       λl+1∕2,j+1 ∕2 4Δ∘ ξ0Δp-----     ξp,l+1∕2,i+1∕2,j
                                     1 - ξ2
                       - --pi+1∕2----------0,j+1λl+1∕2,j+1D (0)           (F.65)
                         λl+1∕2,j+1∕2  4Δ ξ0Δp             ξp,l+1∕2,i+1∕2,j+1
                                   ∘  -------
---                                   1- ξ2
M-(0)[2u]            =   - --pi+1∕2---------0,jλl+1∕2,jD (0)
  p,l+1∕2,i+3∕2,j-1∕2        λl+1∕2,j+1∕2 4Δ ξ0Δp          ξp,l+1∕2,i+1∕2,j
                         ∘ ----2------
                         --1--ξ0,j+1∕2-     (0)
                       -   4Δp Δ ξ0  pi+1D pξ,l+1∕2,i+1,j+1∕2          (F.66)
                         ∘ -----------
--(0)[2u]                    1- ξ20,j+1∕2
M p,l+1∕2,i+1∕2,j-1∕2  =   - ------------pi+1D (0p)ξ,l+1∕2,i+1,j+1∕2
                         ∘ 4Δp-Δ-ξ0---
                           1- ξ2
                       + ------0,j+1∕2piD (0)                    (F.67)
                           4Δp Δ ξ0      pξ,l+1∕2,i,j+1∕2
                                  ∘ ----2--
--(0)[2u]                --pi+1∕2-----1--ξ0,j- l+1∕2,j  (0)
M p,l+1∕2,i- 1∕2,j- 1∕2  =   λl+1∕2,j+1∕2 4Δ ξ0Δp  λ     D ξp,l+1∕2,i+1∕2,j
                        ∘  ----2------
                        ---1--ξ0,j+1∕2    (0)
                       +   4Δp Δξ0   piDpξ,l+1∕2,i,j+1∕2             (F.68)

As expected, both methods merge for uniform momentum and pitch-angle grids, and matrix coefficients for cross-derivatives are similar.