$n+1,\dots, N$.\
Concerning the approximation of the Fokker-Planck equation, following [@acd], we note that the weak formulation of previous equation involves the identity $$-\int_{{\mathbb{T}}^2} {{\rm{div}}}\left(m   H_p(x,D u)\right) w\, dx= \int_{{\mathbb{T}}^2} m H_p(x,D u) \cdot D w\, dx,$$ for any test function $w$.\
Using discrete integration by parts, it can be approximated by $$\label{integralterm}
-h^2 \sum_{i,j} {\mathcal{B}}_{i,j}(U,M) W_{i,j}= h^2\sum_{i,j}    M_{i,j} g_q (x_{i,j},[D_h U]_{i,j}) \cdot [D_h W]_{i,j},$$ where $${\mathcal{B}}_{i,j}(U,M) = 
{\small \frac{1}{h} \left(
\begin{aligned}
&  \left( \begin{aligned} & M_{i,j} \frac{\partial g}{\partial q_1}(x_{i,j}, [D_hU]_{i,j}) - M_{i-1,j} \frac{\partial g}{\partial q_1}(x_{i-1,j}, [D_hU]_{i-1,j}) \\
& + M_{i+1,j} \frac{\partial g}{\partial q_2}(x_{i+1,j}, [D_hU]_{i+1,j}) - M_{i,j} \frac{\partial g}{\partial q_2}(x_{i,j}, [D_hU]_{i,j})  \end{aligned}\right) \\
& + \left( \begin{aligned} & M_{i,j} \frac{\partial g}{\partial q_3}(x_{i,j}, [D_hU]_{i,j}) - M_{i,j-1} \frac{\partial g}{\partial q_3}(x_{1,j-1}, [D_hU]_{1,j-1}) \\
& + M_{i,j+1} \frac{\partial g}{\partial q_4}(x_{i,j+1}, [D_hU]_{i,j+1}) - M_{i,j} \frac{\partial g}{\partial q_4}(x_{i,j}, [D_hU]_{i,j})  \end{aligned}\right)
\end{aligned}\right).}$$ This yields the following scheme: $$\label{discrete_FPequation}
\left\{
\begin{array}{ll}
D_{\Delta t}^{\alpha} M^{n+1}_{i,j}-(\Delta_h M^{n+1})_{i,j} - {\mathcal{B}}_{i,j}(U^{n},M^{n+1})= 0, \quad &n=0,\dots, N-1,\, 0\le i,j\le N_h,\\[6pt]
M_{i,j}^{0}=  \frac 1 {h^2} \int_{|x-x_{i,j}|_{\infty}\le h/2} m_0(x) dx&0\le i,j\le N_h.\\
\end{array}
\right.$$ for a given vector $U=  (U^n_{i,j})_{ 1\le n\le N}$.\
Collecting and, we end up with the following scheme for $$\label{eq:scheme}
\left\{
 \begin{array}{ll}
\bar D_{\Delta t}^{\alpha} U_{i,j}^{n} - (\Delta_h U^{n})_{i,j} + g(x_{i,j},[D_hU^{n}]_{i,j}) = 