---
abstract: 'In this paper, we study the numerical approximation of a system of PDEs with fractional time derivatives. This system is derived from an optimal control problem for a time-fractional Fokker-Planck equation with time dependent drift by convex duality argument. The system is composed by a time-fractional backward Hamilton-Jacobi-Bellman and a forward Fokker-Planck equation and can be used to describe the evolution of probability density of particles trapped in anomalous diffusion regimes. We approximate Caputo derivatives in the system by means of L1 schemes and the Hamiltonian by finite differences. The scheme for the Fokker-Planck equation is constructed such that the duality structure of the PDE system is preserved on the discrete level. We prove the well posedness of the scheme and the convergence to the solution of the continuous problem.'
author:
- Fabio Camilli
- Serikbolsyn Duisembay
- Qing Tang
title: 'Approximation of an optimal control problem for the time-fractional Fokker-Planck equation'
---

[**AMS subject classification**]{}:

:   65N06, 91A13, 35R11 .

[**Keywords**]{}:

:   Mean Field Games, time-fractional Fokker-Planck equation, Caputo derivative, finite differences, convergence.

Introduction
============

In this paper, we study the numerical approximation of the system $$\label{eq:frac_mfg}
\begin{cases}
\partial_{[t,T)}^\alpha u - \Delta u + H(x,Du) = f[m](x) &(x,t)\in{\mathbb{T}}^2\times (0,T),\\[4pt]
\partial_{(0,t]}^\alpha m - \Delta m - {{\rm{div}}}(mH_p(x,Du)) = 0&(x,t)\in{\mathbb{T}}^2\times (0,T),\\[4pt]
u(x,T)=u_T(x),\,m(x,0)=m_0(x), \\[4pt]
\int_{{\mathbb{T}}^2}m_0(x)dx=1,
\end{cases}$$ where ${\mathbb{T}}^2$ denotes the unit torus in ${\mathbb{R}}^2$ and $\partial_{(0,t]}^\alpha\cdot$, $\partial_{[t,T)}^\alpha\cdot$ are the forward and backward time-fractional Caputo derivatives defined, for $\alpha\in (0,1)$, by $$\begin{aligned}
&\partial_{(0,t]}^\alpha m(x,t) = \frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{\partial_s m(x,s)}{(t-s)^{\alpha}}ds, \label{eq:caputo_forward}\\
&\partial_{[t,T)}^\alpha u(x,t) = -\frac{1}{\Gamma(1-\alpha)}\int_{t}^{T}\frac{\partial_su(x,s)}{(s-t)^{\alpha}}ds.\label{eq:caputo_backward}\end{aligned}$$ The system is composed by time-fractional backward Hamilton-Jacobi-Bellman and forward Fokker-Planck equations with the corresponding final and initial conditions. If the coupling cost $f[m]$ admits a primitive $F(x,m)$, can be derived via convex duality from an optimal control problem of $$\inf_{v,m}\int_0^t\int_{{\mathbb{T}}^2} (mL(x,-v)+F(x,m))dx\,ds+\int_{{\mathbb{T}}^2}(I^{1-\alpha}_{(0,t]}m)(T)dx,$$ subject to the constraint given by the time-fractional Fokker-Planck equation $$\label{frac_fp}
\begin{cases}
\partial_{(0,t]}^\alpha m - \Delta m - {{\rm{div}}}(mv) = 0,&(x,t)\in{\mathbb{T}}^2\times (0,T),\\[4pt]
m(x,0)=m_0(x),\\[4pt]
\int_{{\mathbb{T}}^2}m_0(x)dx=1.
\end{cases}$$ Here the forward Riemann-Liouville fractional integral is defined as $$I^{1-\alpha}_{(0,t]}m(t)=\frac{1}{\Gamma(1-\alpha)}\int_0^t\frac{m(s)}{(t-s)^\alpha}ds.$$ System is derived in [@q], where the Fokker-Planck equation is considered in its weak formulation via a general fractional Sobolev space framework. The wellposedness of the coupled system is obtained using compactness and Schauder fixed-point argument. Time fractional advection-diffusion equations have been extensively used for studying evolution of probability of particles governed by subdiffusion behavior. This refers to the trajectory of a single particle is modeled by non Markovian stochastic process where the mean square displacement is no longer linear [@abmw; @bg; @mk; @ms; @zaslavsky2002]. The nonlocal structure of fractional time derivative makes it very useful for describing memory effects and deriving the power law decay structure of energy in a diffusion-transport PDE system [@du; @jlz; @Jin2019; @tt].\
Time fractional Hamilton-Jacobi-Bellman equation has been considered from a probabilistic point of view by Kolokoltsov et al.[@Kolokoltsov2014]. The theory of viscosity solution of time-fractional Hamilton-Jacobi-Bellman equation has been constructed by Topp and Yangari [@Topp2017], Giga, Namba et al. [@glm; @Namba2018]. Recently Camilli and Goffi studied weak solutions to time-fractional Hamilton-Jacobi-Bellman equation in a Sobolev space [@cg]. The fractional Fokker-Planck equation can be more suitably considered in the variational formulation and the weak solution is in the sense of distributions. The natural functional structure would be that of vector valued fractional Sobolev space and wellposedness is usually established by a priori estimates based on energy dissipation. This approach has been adopted in a large number of works, [@du; @du2018invitation; @Gorenflo2015; @Jin_2015; @Jin2019; @Shen_2019; @q; @tt] for examples.\
From a numerical point view, the essential step to study the fractional PDE is to approximate the time-fractional derivative operator. We use a classical scheme called L1 approximation. This scheme is naturally derived from the approximation of the fractional integral as a Riemann sum and long known to be consistent (for an concise illustration we refer to section 3 of [@tt]). Moreover, L1 approximation scheme stands out by being able to preserve at the discrete level certain desirable features of the origin PDEs, such as maximum principle [@glm] and energy stability [@du; @Shen_2019; @tt]. Therefore, approximations can be constructed to converge to suitable notions of weak solutions, which may be viscosity solutions to Hamilton-Jacobi-Bellman equations [@glm], or distributional solutions in Sobolev spaces to diffusion (phase field) equations [@du; @Shen_2019; @tt].\
If fractional derivatives are simply replaced by first order derivatives, system becomes the classical Mean Field Games system introduced in [@ll] to characterize the Nash equilibria for differential games with a large number of agents. Achdou et al. [@acd; @acc] introduced an finite difference type discretization of mean field games system where the scheme for the Fokker-Planck equation is deduced via a duality argument, reproducing at a discrete level a crucial property of the continuous system. This type of schemes has been further developed in [@afg; @bks; @cs].\
In this paper, we present a finite differences approximation of by combining the classical approximation scheme for mean field games system developed in [@acd; @acc] with L1 schemes for the fractional derivatives. The duality structure of the time-fractional PDE system and the energy stable structure of the advection diffusion equation are shown to be preserved. Under appropriate monotonicity and convexity assumptions, we prove existence, uniqueness of the solution to the discrete system and convergence to the solution of the original PDE system. Some numerical tests are performed to illustrate the different behavior of the time-fractional case with respect to the standard one.\
The paper is organized as follows. In Section \[sec:scheme\], we introduce the finite differences scheme and prove some preliminary results. Section \[sec:well\_posed\] is devoted to show the well posedness of the discrete problem. In Section \[sec:convergence\], we prove the convergence of the scheme and in Section \[sec:numerics\] we carry out some numerical simulations. Finally we discuss some possible directions for future works.

The numerical scheme {#sec:scheme}
====================

In this section we introduce a numerical scheme for and we prove some preliminary results to show its well-posedness.\
For simplicity of notations, we assume that the dimension $d$ is equal to $2$, since the extension for general $d$ will be clear from this special case. Let ${\mathbb{T}}^2_h$ be a uniform grid on the torus with step $h$, (this supposes that $N_h = 1/{h}$ is an integer), and denote by $x_{i,j}$ a generic point in ${\mathbb{T}}^2_h$ (an anisotropic mesh with steps $h_1$ and $h_2$ is possible too and we have taken $h_1=h_2$ only for simplicity). The values $U^n_{i,j}$ and $M^n_{i,j}$ denotes the numerical approximation of the function $u$ and $m$ at $(x_{i,j},t_n)=(ih,jh, n{{\Delta t}})$, $0\le i,j\le N_h$, $n=0,\dots,N$ (assuming that $N=T/{{\Delta t}}$ is an integer). We also denote by $U^n$ and $M^n$ the grid function taking the values $U^n_{i,j}$ and $M^n_{i,j}$ at $x_{i,j}\in {\mathbb{T}}^2_h$. For a generic pair of grid functions $U$, $V$ on ${\mathbb{T}}^2_h$, we define the scalar product\
$$(U,V)_2=h^2\sum_{  i,j } U_{i,j} V_{i,j},$$ the norm\
$$\|U\|_2^2=(U,U)_2$$ and the norm $$\|U\|_\infty=\max_{0\le i,j <N_h}|U_{ij}|.$$ We consider the compact and convex set $${\mathcal{K}}_h=\{M= (M_{i,j})_{ 0\le i,j <N_h}:   (M,\overline 1)_2=1;\quad   M_{i,j}\ge 0 \}\,,$$ where $\overline 1$ denotes the $N_h^2$-tuple with all components equal to $1$. Note that ${\mathcal{K}}_h$ can be viewed as the set of the discrete probability measures on ${\mathbb{T}}^2_h$.\
We first recall the L1 numerical method for the discretization of the forward and backward Caputo time-fractional derivatives (see [@jlz; @lx]). The approximations of and at $(x_{i,j},t_n)$ are given by $$\begin{aligned}
D_{\Delta t}^{\alpha} M_{i,j}^{n} = \frac{1}{\rho_\alpha} \left( M_{i,j}^{n} - \sum_{k=0}^{n-1} c_k^{n} M_{i,j}^{k} \right),\label{eq:discr_frac_deriv_forw}\\
\bar D_{\Delta t}^{\alpha} U_{i,j}^{n} = \frac{1}{\rho_\alpha}\left( U_{i,j}^{n} - \sum_{k=n+1}^{N} {\overline c}_{n}^k U_{i,j}^k\right) ,\label{eq:discr_frac_deriv_back}\end{aligned}$$ where $$\rho_\alpha = \Gamma(2-\alpha)\Delta t^{\alpha}$$ and the forward and backward coefficients are defined by $$c_k^n = 
\begin{cases}
n^{1-\alpha}-(n-1)^{1-\alpha}, \quad k=0,\\
2(n-k)^{1-\alpha} - (n+1-k)^{1-\alpha} - (n-1-k)^{1-\alpha}, \quad k=1,\dots, n-1,
\end{cases}$$ and $${\overline c}_n^k = 
\begin{cases}
2(k-n)^{1-\alpha} - (k+1-n)^{1-\alpha} - (k-1-n)^{1-\alpha}, \quad k=n+1,\dots, N-1, \\
(N-n)^{1-\alpha}-(N-1-n)^{1-\alpha}, \quad k=N.
\end{cases}$$ In the following lemma, we summarize some properties of the coefficients $c_k^n$, ${\overline c}_n^k$.

\[lem:c\_properties\] We have

1.  $c_k^{n} >0$ for $0\le k \leq n-1$ and ${\overline c}_n^m >0$ for $ n+1\le m\le N$.

2.  $\bar c_k^n = c_{k+1}^{n+1}$ for $1\le n\le N-1$ and $0\le k\le n-1$.

3.  $\sum_{k=0}^{n-1} c_k^{n}=1$ and $\sum_{k=n+1}^{N} {\overline c}_n^k=1.$

By simple calculation, it can be shown that an equivalent form of the L1 scheme of the discrete forward Caputo derivative is $$\label{discreteCaputo2}
D_{\Delta t}^{\alpha} U_{i,j}^{n} =\sum_{k=0}^nb_k\frac{U_{i,j}^{n-k}-U_{i,j}^{n-1-k}}{\Delta t},$$ where $$b_k=\frac{\Delta t^{1-\alpha}}{\Gamma(2-\alpha)}[(k+1)^{1-\alpha}-k^{1-\alpha}].$$ When written in this form and used for fractional phase field equations it has been shown the L1 scheme of Caputo derivative can preserve on a discrete level maximum principle, energy stability and mass conservation [@tt].\
For the discretization in space, we rely on the approach in [@acd]. We introduce the finite differences operators $$\label{eq:finite_diff}
(D_1^+ U )_{i,j} = \frac{ U_{i+1,j}-U_{i,j}   } {h} \quad \hbox{and }\quad  (D_2^+ U )_{i,j} = \frac{ U_{i,j+1}-U_{i,j}   } {h},$$ and define $$[D_h U]_{i,j} =\left((D_1^+ U )_{i,j} , (D_1^+ U )_{i-1,j}, (D_2^+ U )_{i,j}, (D_2^+ U )_{i,j-1}\right) ^T.\label{eq:discrete_gradient}$$ For the discretization of the Laplace operator, we consider a standard five points approximation $$(\Delta_h U)_{i,j} = -\frac{1}{h^2} (4 U_{i,j} - U_{i+1,j} - U_{i-1,j} - U_{i,j+1} - U_{i,j-1}).$$ In order to approximate the Hamiltonian $H$ in equation , we consider a numerical Hamiltonian $g: {\mathbb{T}}^2 \times {\mathbb{R}}^4\to {\mathbb{R}}$, $(x,q_1,q_2,q_3,q_4)\mapsto g\left(x,q_1,q_2,q_3,q_4\right)$ satisfying the following conditions:

-   $g$ is nonincreasing with respect to $q_1$ and $q_3$, and nondecreasing with respect to $q_2$ and $q_4$.

-   $g$ is consistent with the Hamiltonian $H$, i.e. $$g(x,q_1,q_1,q_2,q_2)=H(x,q), \quad \forall x\in {\mathbb{T}}^2,   q=(q_1,q_2)\in {\mathbb{R}}^2.$$

-   $g$ is of class $C^1$.

-   the function $(q_1,q_2,q_3,q_4) \mapsto g(x,q_1,q_2,q_3,q_4)$ is convex.

For the discretization of the coupling term $f[m]$ in the Hamilton-Jacobi-Bellman equation, denoted by $m_h$ the piecewise constant function with value $M_{i,j}$ in the square $|x-x_{i,j}|_\infty \le h/2$, we assume that $f[m_h]$ can be calculated in practice. Then, the approximation of $f[m](x_{i,j})$ is given by $(f_h[M])_{i,j} = f[m_h](x_{i,j})$. In particular, if $f$ is a local operator, i.e. $f[m](x)=f(m(x))$, then we set $f_h[M]_{j,k}=f(M_{j,k})$.\
We assume that:

-   $f_h$ is continuous and maps ${\mathcal{K}}_h$ on a bounded set of grid functions.

-   $f_h$ is monotone, i.e. for $M,\bar M\in {\mathcal{K}}_h$ $$\left( f_h[M]-f_h[\bar M],M-\bar  M \right)_2 \le 0 \Rightarrow  M=  \bar M.$$

-   There exists $C$ independent of $h$ such that for all grid functions $M\in {\mathcal{K}}_h$ $$\begin{aligned}
        & \| f_h[M] \|_{\infty} \le C\\
        & | f_h[M]_{i,j} - f_h[M]_{\bar i,\bar j} | \le C |x_{i,j}-x_{\bar i,\bar j}|\qquad \forall 0\le i,j,\bar i,\bar j \le N_h.
        \end{aligned}$$

-   For any $m\in {\mathcal{K}}=\{m\in L^1({\mathbb{T}}^2): m\ge 0, \int_{{\mathbb{T}}^2}m(x)dx=1\}$, defined $$M_h=\frac 1 {h^2} \int_{|x-x_{i,j}|_{\infty}\le h/2} m(x)dx,$$ then $$\lim_{h\to 0}\,\sup_{m\in{\mathcal{K}}}\|f[m]-f_h[M_h]\|_\infty =0.$$

We approximate the Hamilton-Jacobi-Bellman equation in by means of the scheme $$\label{discrete_HJequation}
\left\{
\begin{array}{ll}
\bar D_{\Delta t}^{\alpha} U^n_{i,j} - (\Delta_h U^{n})_{i,j} + g(x_{i,j},[D_h U^{n}]_{i,j}) = f_h[M^{n+1}]_{i,j}, \quad& n=0,\dots, N-1,\, 0\le i,j\le N_h,\\[6pt]
U^N_{i,j}=u_T(x_{i,j})& 0\le i,j\le N_h,
\end{array}
\right.$$ for a given $M= (M^n)_{ 0\le n\le N-1}$. Note that the previous scheme is backward, implicit and non local in the time variable, since the solution at step $n$ depends on the values of the solution computed at previous steps $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}) = f_h[M^{n+1}]_{i,j} ,\\[6pt]
D_{\Delta t}^\alpha M_{i,j}^{n+1} - (\Delta_h M^{n+1})_{i,j} - {\mathcal{B}}_{i,j}(U^{n},M^{n+1}) = 0 ,\\[6pt]
U^{N_T}_{i,j} = u_T(x_{i,j}), \,   M_{i,j}^{0}=  \frac 1 {h^2} \int_{|x-x_{i,j}|_{\infty}\le h/2} m_0(x) dx,
\end{array}
\right.$$ for $n=0,\dots,N-1$ and $0\le i,j\le N_h$.\
We prove some results for the problems and which will be used in the following.

\[le:existence\] Let $V$ be a grid function on ${\mathbb{T}}_h^2$ and $\delta$ be a positive parameter. Then, there exists a unique grid function $U$ such that $$\label{eq:stationary}
\delta U_{i,j} - (\Delta_hU)_{i,j} + g(x_{i,j},[D_hU]_{i,j}) = V_{i,j}.$$

For the proof of the previous result, we refer to [@acd Lemma 1].

\[prop:ex\_HJ\] There exists a function $U=(U^n)_{n=0,\dots,N}$ which satisfies . Moreover, if $U= (U^n)_{ 0\le n\le N }$, $V= (U^n)_{ 0\le n\le N }$ are such that $U^N\le V^N$ and $$\begin{aligned}
\bar D_{\Delta t}^{\alpha} U_{i,j}^{n} - (\Delta_h U^{n})_{i,j} + g(x_{i,j},[D_hU^{n}]_{i,j}) \le f_h[M^{n+1}]_{i,j},\quad n=0,\dots,N-1,\\
\bar D_{\Delta t}^{\alpha} V_{i,j}^{n} - (\Delta_h V^{n})_{i,j} + g(x_{i,j},[D_hV^{n}]_{i,j}) \ge  f_h[M^{n+1}]_{i,j},\quad n=0,\dots,N-1,\end{aligned}$$ then $U^n\le V^n$ for any $n=0,\dots,N$.

The existence and uniqueness of $U^{n}$, $n=0, \dots, N-1$ is obtained by induction: at each step, we use Lemma \[le:existence\] with $$\delta = \frac{1}{\rho_\alpha},\qquad V_{i,j} =  \frac{1}{\rho_\alpha} \sum_{k=n+1}^{N} \bar c_{n}^k U_{i,j}^k + f_h[M^{n+1}]_{i,j}.$$ The second part of the statement follows immediately by a comparison principle for , which relies on the monotonicity of $g$ and the positivity of the coefficients ${\overline c}_n^k$.

To prove that a solution of is uniformly Lipschitz continuous, we need the following fractional version of the discrete Gronwall Lemma (see [@lwz Lemma 3.1]).

Suppose that the nonnegative sequences $\{{\omega}^n,g^n\}_{n\in{\mathbb{N}}}$ satisfy $D_{\Delta t}^\alpha{\omega}^1\le {\lambda}_1 {\omega}^1+{\lambda}_2 {\omega}^0+g^1$ and $$D_{\Delta t}^\alpha {\omega}^n\le {\lambda}_1{\omega}^n+{\lambda}_2{\omega}^{n-1}+{\lambda}_3{\omega}^{n-2}+g^n,\quad n\ge 2$$ where ${\lambda}_1,{\lambda}_2,{\lambda}_3$ are given nonnegative constants independent of ${{\Delta t}}$. Then, there exists $\delta>0$ such for ${{\Delta t}}<\delta$, $$\label{eq:discrete_Gronwall}
{\omega}^n\le \left({\omega}^0+\frac{(n{{\Delta t}})^\alpha}{\Gamma(1+\alpha)}\max_{1\le j \le n}g^j\right)E_\alpha(2\lambda (n{{\Delta t}})^\alpha ),\qquad 1\le n\le N,$$ where $E_\alpha(z)=\sum_{k=0}^\infty z^k/\Gamma(1+k\alpha)$ is the Mittag-Leffler function and ${\lambda}={\lambda}_1+{\lambda}_2/(2-2^{1-\alpha})+{\lambda}_3/(2\cdot2^{1-\alpha}-1^{1-\alpha}-3^{1-\alpha})$.

\[prop:lip\_U\] Let $U=(U^n)_{n=0,\dots,N}$ be the solution of . Then, for ${{\Delta t}}$ sufficiently small, there exists a positive constant $L$ such that $$\label{eq:lip_est}
    \|D_hU^n\|_\infty \le L, \qquad  n=0,\dots,N.$$

Rewrite as the stationary equation $$U_{i,j}^{n} - \rho_\alpha (\Delta_h U^{n})_{i,j} + \rho_\alpha g(x_{i,j},[D_hU^{n}]_{i,j}) - \rho_\alpha f_h[M^{n+1}]_{i,j}= \sum_{k=n+1}^{N} {\overline c}_{n}^k U_{i,j}^k$$ and let $\Psi$ be the map $U_{i,j}^n=\Psi( \sum_{k=n+1}^{N} {\overline c}_{n}^k U_{i,j}^k)$ that associates to right hand side of the previous equation the corresponding solution. For $V \in {\mathbb{R}}^{N_h}$ and $\lambda\in{\mathbb{R}}$, we have $\Psi( V+\lambda)=\Psi( V)+\lambda$ (where, on the left side of the previous identity, $\lambda$ is interpreted as the constant function on ${\mathbb{R}}^{N_h}$) and, by standard monotonicity argument, $$\label{eq:lip0}
\|\Psi(V)-\Psi(W)\|_\infty\le  \|V-W\|_\infty.$$ for any $V, W\in {\mathbb{R}}^{N_h}$. For fixed $l,\kappa \in {\mathbb{Z}}^2$, call $\tau U$ the discrete function defined by $$(\tau U)^n_{i,j} = U^n_{l+i,\kappa+j}.$$ We have that $$\bar D_{\Delta t}^{\alpha} (\tau  U)_{i,j}^{n} - (\Delta_h (\tau U)^{n})_{i,j} + g(x_{i,j},[D_h(\tau U)^{n}]_{i,j})= f_h[M^{n+1}]_{i,j}+ E_{i,j}$$ where $$\begin{aligned}
E_{i,j}& =   f_h[M^{n+1}]_{i+l,j+\kappa} -f_h[M^{n+1}]_{i,j}  - g(x_{i+l,j+\kappa},[D_h(\tau U)^{n}]_{i,j}) + g(x_{i,j},[D_h(\tau U)^{n}]_{i,j}).\end{aligned}$$ Hence $ (\tau U)^n=  \Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k (\tau U)^k+\rho_\alpha E\right)$. Moreover, since $$\begin{aligned}
&f_h[M^{n+1}]_{i+l,j+\kappa} -f_h[M^{n+1}]_{i,j} \le C h \sqrt{l^2+\kappa^2},\\
&g(x_{i,j},[D_h(\tau U)^{n}]_{i,j}) - g(x_{i+l,j+\kappa},[D_h(\tau U)^{n}]_{i,j}) \le C(1+\|D_hU^n\|_\infty)h\sqrt{l^2+\kappa^2}\end{aligned}$$ it follows that $$\|E\|_\infty \le C(1+\|D_hU^n\|_\infty)h\sqrt{l^2+\kappa^2}.$$ By , we have $$\begin{aligned}
&\|(\tau U)^n-U^n\|_\infty= \left\|\Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k (\tau U)^k+\rho_\alpha E\right)
-\Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k  U^m\right)\right\|_\infty\\
&\le \left\|\sum_{k=n+1}^{N} {\overline c}_{n}^k( (\tau U)^k - U^k)\right\|_\infty+{\rho_\alpha}\|E\|_\infty.\end{aligned}$$ Hence $$\left\|\frac{(\tau U)^n - U^n}{h \sqrt{l^2+\kappa^2} }\right\|_\infty \le   \sum_{k=n+1}^{N}{\overline c}_{n}^k\left\| \frac{(\tau U)^k - U^k}{h \sqrt{l^2+\kappa^2}}\right\|_\infty + C \rho_\alpha (1+ \|D_hU^{n}\|_\infty).$$ For the arbitrariness of $l,\kappa$, the previous inequality implies that $$\bar D_{\Delta t}^{\alpha}(\|D_hU^n\|_\infty)\le C\|D_hU^n\|_\infty+C$$ and therefore, applying a backward version of the discrete Gronwall’s inequality with ${\omega}^n=\|D_hU^n\|_\infty$, ${\lambda}_1=C$, ${\lambda}_2={\lambda}_3=0$ and $g^n=C$, we obtain for some positive constant $L$.

\[prop:ex\_FP\] For ${{\Delta t}}$ sufficiently small, there exists a unique solution $(M^n)_{n=0}^N$ to . Moreover, $M^n\ge 0$ and $(M^n,\bar 1)_2=1$ for any $n=0,\dots,N$.

Equation is equivalent to $$\label{eq:FP_0} 
    M^{n+1}_{i,j} -\rho_\alpha(\Delta_h M^{n+1})_{i,j} - \rho_\alpha{\mathcal{B}}_{i,j}(U^{n},M^{n+1})  =  \sum_{k=0}^{n} c_k^{n+1} M_{i,j}^{k}.$$ For fixed $n$, is a linear problem in the unknown $M^{n+1}\in ({\mathbb{R}}^{N_h^2})$, which can be rewritten as $$\label{eq:FP_1}
    (I+\rho_\alpha A)M^{n+1}=\sum_{k=0}^{n} c_k^{n+1} M^{k}.$$ where $A$ is the linear operator given by $A\cdot =-  \Delta_h\cdot - {\mathcal{B}}_{i,j}(U^{n},\cdot)$. We observe that, due to the monotonicity assumption (G1), the matrix corresponding to $A$ has positive diagonal entries and non positive off-diagonal entries. Hence, for ${{\Delta t}}$ sufficiently small, the matrix corresponding to the linear operator $I+\rho_\alpha A$ in is a $M$-matrix and therefore non singular. It follows that, for any $n=0,\dots,N-1$, there exists a unique solution to . Moreover, since $M^0\ge 0$ and the inverse of a $M$-matrix is non negative, we easily get by induction that $M^n\ge 0$ for $n=0,\dots,N-1$.\
To show that $(M^n,\bar 1)_2=1$ for any $n=0,\dots,N$, we first recall the following identity (see [@acd eq. (50)]) $$\label{eq:bypart_space}
\begin{aligned}
(AW,Z)_2&=\sum_{i,j}(D^+_1W)_{i,j}(D^+_1Z)_{i,j}+\sum_{i,j}(D^+_2W)_{i,j}(D^+_2Z)_{i,j}\\
&+\sum_{i,j}W_{i,j}[D_h Z]_{i,j}\cdot D_q g(x_{i,j},[D_hU]_{i,j}).
\end{aligned}$$ From the previous identity, it follows that $(AW,\bar 1)_2=0$ for any grid function $W$. Hence, taking the inner product of with the function $\bar 1$, we get that $$(M^{n+1},\bar 1)_2=\sum_{k=0}^{n } c_k^{n+1} (M^{k},\bar 1)_2.$$ Recalling that $M^0\in {\mathcal{K}}_h$ and Lemma \[lem:c\_properties\].(iii), by the previous formula we get by induction that $M^{n} \in {\mathcal{K}}_h$ for any $n=0,\dots,N$.

We can provide an alternative proof to the discrete mass conservation property $(M^n,\bar 1)_2=1$ by the following observation. Since for integrating on the torus ${\mathbb{T}}^2$ the second equation, noting $$\int_{{\mathbb{T}}^2}\Delta m-{{\rm{div}}}(mH_p(x,Du))dx=0,$$ so that $$\begin{aligned}
0&=\int_{{\mathbb{T}}^2}\partial_{(0,t]}^{\alpha}m(x,t)dx=\frac{1}{\Gamma(1-\alpha)}\int_{{\mathbb{T}}^2}\int_0^t\frac{1}{(t-s)^{\alpha}}\frac{\partial m(x,s)}{\partial s}dsdx\\
&=\frac{1}{\Gamma(1-\alpha)}\int_0^t\frac{1}{(t-s)^{\alpha}}\big{(}\int_{{\mathbb{T}}^2}\frac{\partial m(x,s)}{\partial s}dx\big{)}ds,\,\,\,\,\,\forall t\in [0,T].\end{aligned}$$ Then $$\frac{d}{dt}\int_{{\mathbb{T}}^2}m(x,t)dx=\int_{{\mathbb{T}}^2}\frac{\partial}{\partial t}m(x,t)dx=0,\,\,\,\,\forall t\in [0,T].$$ Therefore we have $$\int_{{\mathbb{T}}^2}m(x,t)dx=\int_{{\mathbb{T}}^2}m(x,0)dx=1,\,\,\,\,\forall t\in [0,T].$$ Analogously, for the discrete case we can obtain $$\sum_{i,j}(\Delta_h M^{n+1})_{i,j} +{\mathcal{B}}_{i,j}(U^{n},M^{n+1})=0$$ and $$\begin{aligned}
0=h^2\sum_{i,j}D^{\alpha}_{\Delta t}M_{i,j}^{n+1}=D^{\alpha}_{\Delta t}(M^{n+1},\bar{1})_2,\,\,\,0\leq n\leq N-1.\end{aligned}$$ By using the equivalent form of L1 scheme for forward Caputo derivative we have $$\sum_{k=0}^{n}b_k\frac{(M^{n+1-k},\bar{1})_2-(M^{n-k},\bar{1})_2}{\Delta t}=0,\,\,\,0\leq n\leq N-1.$$ It follows easily by mathematical induction that $$(M^{n},\bar{1})_2=(M^{0},\bar{1})_2=1,\,\,\,1\leq n\leq N.$$

The following lemma provides an integration by parts for , which is the discrete counterpart of fractional integration by parts formula with Caputo derivative (see [@q Lemma 2.2]) $$\int_{0}^T \partial_{(0,t]}^{\alpha} \phi(t)\kappa(t)dt + \phi(0)(I_{[t,T)}^{1-\alpha}\kappa)(0) = \int_{0}^T \phi(t) \partial_{[t,T)}^{\alpha} \kappa(t)dt + \kappa(T)(I_{(0,t]}^{1-\alpha}\phi)(T),$$ where $\phi(t), \kappa(t)\in C^1[0,T]$. \[le:by\_parts\_n\_n\_plus\_1\] Given $U=\{U^n\}_{n=0}^{N}$, $M=\{M^n\}_{n=0}^{N}$, then $$\begin{aligned}
&\sum_{n=0}^{N-1} \Big(\bar D_{\Delta t}^{\alpha} U^{n}, M^{n+1}\Big)_2 + \frac{1}{{\rho_\alpha}}\sum_{n=0}^{N-1} \bar c_n^N \left(U^N,M^{n+1}\right)_2  \\
&= \sum_{n=0}^{N-1} (D_{\Delta t}^{\alpha} M^{n+1},U^{n})_2+ \frac{1}{{\rho_\alpha}}\sum_{n=0}^{N-1} c_0^{n+1}\left(M^0,U^n\right)_2.\end{aligned}$$ For the proof, see the Appendix. We also recall the following result [@glm Lemma 2.2], which will be used for barrier arguments.

\[le:barrier\_argument\] For ${{\Delta t}}>0$, the following inequality holds $$\label{eq:barrier1}
    \bar D_{\Delta t}^\alpha ((N-n){{\Delta t}})^\alpha\ge \frac{\alpha(1-\alpha)}{\Gamma (2-\alpha)},\qquad n=0,\dots,N-1.$$

Well posedness of the discrete problem {#sec:well_posed}
======================================

In this section, we study the well posedness of the discrete problem .

Uniqueness
----------

The uniqueness argument is classical in Mean Field Games theory and relies on the monotonicity of the coupling cost $f$ and convexity of the discrete Hamiltonian $g$.

\[thm:uniqueness\] The problem has at most one solution.

We suppose by contradiction that there are two solutions $(U,M)$ and $(\tilde{U},\tilde{M})$ such that $U^N = \tilde{U}^N$ and $M^{0} = \tilde{M}^{0}$. Then we have $$\begin{aligned}
\bar D_{\Delta t}^{\alpha} (U^{n}& -\tilde{U}^{n})_{i,j} - (\Delta_h (U^{n}-\tilde{U}^{n}))_{i,j} \\
& + g(x_{i,j},[D_hU^{n}]_{i,j}) - g(x_{i,j},[D_hU^{n}]_{i,j}) -  (f_h[M^{n+1}] - f_h[\tilde{M}^{n+1}] )_{i,j}=0.
\end{aligned}$$ Multiplying by $Z_{i,j}$ and summing over all $i,j$ leads to $$\label{eq:u_diff}
\begin{aligned}
(\bar D_{\Delta t}^{\alpha} & (U^{n} -\tilde{U}^{n}),Z)_2 - ((\Delta_h (U^{n}-\tilde{U}^{n})),Z)_2 \\
 & + \sum_{i,j} (g(x_{i,j},[D_hU^{n}]_{i,j}) - g(x_{i,j},[D_h\tilde{U}^{n}]_{i,j}))Z_{i,j} - (( f_h[M^{n+1}] - f_h[\tilde{M}^{n+1}]),Z )_2=0.
\end{aligned}$$ Multiplying the discrete equation by $W_{i,j}$ and summing over all $i,j$ we obtain $$(D_{\Delta t}^\alpha M^{n+1},W)_2 - (M^{n+1}, \Delta_h W)_2 + \sum_{i,j} M_{i,j}^{n+1} [D_hW]_{i,j} \cdot \nabla_q g  (x_{i,j},[D_hU^{n}]) = 0.$$ From this and the similar equation satisfied by $\tilde{M}^{n+1}$, we obtain $$\label{eq:m_diff}
\begin{aligned}
0 & = (D_{\Delta t}^\alpha (M^{n+1}  -\tilde{M}^{n+1}),W)_2 - ((M^{n+1}-\tilde{M}^{n+1}), \Delta_h W)_2 \\
& + \sum_{i,j} M_{i,j}^{n+1} [D_hW]_{i,j} \cdot \nabla_q g ( x_{i,j},[D_hU^{n}])  - \sum_{i,j} \tilde{M}_{i,j}^{n+1} [D_hW]_{i,j} \cdot \nabla_q g  (x_{i,j},[D_h\tilde{U}^{n}] ).
\end{aligned}$$

Let $Z = M^{n+1} - \tilde{M}^{n+1}$ in and $W = U^{n} - \tilde{U}^{n}$ in . Subtracting the two equations so obtained, summing over $n=0,1,\dots,N-1$ and recalling that, by Lemma \[le:by\_parts\_n\_n\_plus\_1\], we have $$\sum_{n=0}^{N-1} (\bar D_{\Delta t}^{\alpha} (U^{n} -\tilde{U}^{n}),M^{n+1} - \tilde{M}^{n+1})_2 =  \sum_{n=0}^{N-1} (D_{\Delta t}^\alpha (M^{n+1}  -\tilde{M}^{n+1}),U^{n} - \tilde{U}^{n})_2,$$ we finally get $$\label{eq:un_nonneg_terms}
\begin{aligned}
0 & = \sum_{n=0}^{N-1}   \bigg(f_h[M^{n+1}] - f_h[\tilde{M}^{n+1}],M^{n+1} - \tilde{M}^{n+1}\bigg)_2 \\
& +   \sum_{n=0}^{N-1}  \sum_{i,j} M_{i,j}^{n+1}\bigg( g(x_{i,j},[D_h\tilde{U}^{n}]_{i,j}) - g(x_{i,j},[D_hU^{n}]_{i,j})  \\
& -[D_h(\tilde{U}^{n}-U^{n})]_{i,j} \cdot \nabla_q g ( x_{i,j},[D_hU^{n}])   \bigg)\\
& + \sum_{n=0}^{N-1}  \sum_{i,j} \tilde{M}_{i,j}^{n+1} \bigg(g(x_{i,j},[D_hU^{n}]_{i,j}) - g(x_{i,j},[D_h\tilde{U}^{n}]_{i,j}) \\
& -[D_h(U^{n}-\tilde{U}^{n})]_{i,j} \cdot \nabla_q g ( x_{i,j},[D_h\tilde{U}^{n}] )  \bigg).
\end{aligned}$$ By the monotonicity of $f_h[\cdot]$ and assumption (G4), all three terms on RHS of are non negative. Hence, we obtain that $M^{n} = \tilde M^{n}$ for all $0 \le n \le N$. By Prop. \[prop:ex\_HJ\], it also follows that $U^{n} = \tilde U^{n}$ for all $0 \le n \le N$.

Existence
---------

To show the existence of a solution for the scheme , we rely on a fixed point argument.

For ${{\Delta t}}$ sufficiently small, has a solution.

First observe that, by the assumption on $m_0$, we have that $M^0\in{\mathcal{K}}_h$. We denote $$\label{eq:def_K_N}
 {\mathcal{K}}_h^N=\{M=(M^n)_{0\le n\le N} :\, M^0= \frac 1 {h^2} \int_{|x-x_{i,j}|_{\infty}\le h/2} m_0(x) dx,\,  M^n\in{\mathcal{K}}_h,\,\forall n=0,\dots,N\}$$ and we define a map $\Phi$ in the following way: Let $\Psi_1:{\mathcal{K}}^N_h\to ({\mathbb{R}}^{N_h^{2}})^N$ be the map which associates to $M=(M^n)_{0\le n\le N}$ the solution $U=(U^n)_{0\le n\le N}$ of and let $\Psi_2: ({\mathbb{R}}^{N_h^2})^N\to ({\mathbb{R}}^{N_h^2})^N$ the map that to $U=(U^n)_{0\le n\le N}$ associates the solution of . Then we show that the map $\Phi=\Psi_2\circ \Psi_1$, is well defined on ${\mathcal{K}}^N_h$ and maps this set into itself.\
We show that the map $\Psi_1$ is bounded by means of a barrier argument. By , the function $$V^n=\|U^N\|_\infty + \frac{M_0\Gamma (2-\alpha)}{\alpha(1-\alpha)}((N-n){{\Delta t}})^\alpha,$$ where $M_0=\sup_{i,j}|H(x_{i,j},0)-(f_h(M^{n+1}))_{i,j}|$, satisfies $V^N\ge U^N$ and $$\bar D_{\Delta t}^{\alpha} V_{i,j}^{n} - (\Delta_h V^{n})_{i,j} + g(x_{i,j},[D_hV^{n}]_{i,j}) \ge  f_h[M^{n+1}]_{i,j},\quad n=0,\dots,N-1.$$ By Proposition \[prop:ex\_HJ\], we get $U^n\le V^n$, $n=0,\dots,N$. In a similar way, we prove that $U^n\ge- V^n$, $n=0,\dots,N$, and therefore $$\|U^n\|_\infty \le \|U^N\|_\infty +CT^\alpha.$$ Hence $\Psi_1$ maps ${\mathcal{K}}_h^N$ in a bounded subset of $({\mathbb{R}}^{N_h^2})^N$.\
To prove that the map $\Psi_1$ is continuous, consider a sequence $M_{\kappa}=(M_{\kappa}^n)_{0\le n\le N}\in {\mathcal{K}}_h^N$ which tends to $M=(M^n)_{0\le n\le N}\in {\mathcal{K}}_h^N$ for $\kappa \to \infty$ and denote by $U_{\kappa}$, $U \in ({\mathbb{R}}^{N_h^2})^N$ the corresponding solutions of . Set $$E_{\kappa} = \sup_{ n=0,\dots,N-1}  \left\| f_h[M_{\kappa}^{n+1}]  -f_h[M^{n+1}] \right\|_\infty$$ and observe that, by assumption (F1) on $f_h$, $E_{\kappa}\to 0$ for $k\to \infty$. Set $V^n_{\kappa}=U^n_{\kappa}- \frac{\Gamma(2-\alpha)}{\alpha(1-\alpha)}((N-n){{\Delta t}})^\alpha E_{\kappa}$, $n=0,\dots,N$, then $$\begin{aligned}
\bar D_{\Delta t}^{\alpha} V^n_{i,j}  - (\Delta_h V^{n})_{i,j} + g(x_{i,j},[D_h V^{n}]_{i,j})\le  f_h[M_k^{n+1}]_{i,j}-E_{\kappa}= \\
f_h[M^{n+1}]_{i,j}+f_h[M_k^{n+1}]_{i,j}-f_h[M^{n+1}]_{i,j}-  E_{\kappa}
\le f_h[M^{n+1}]_{i,j}.\end{aligned}$$ Moreover $V^N\le U^N$ and therefore, by Prop. \[prop:ex\_HJ\], $V^n_{\kappa}\le U^n$ for any $n=0,\dots,N$. Using a similar argument for $W^n_{\kappa}=U^n_{\kappa}+\frac{\Gamma(2-\alpha)}{\alpha(1-\alpha)}((N-n){{\Delta t}})^\alpha E_{\kappa}$, $n=0,\dots,N$, we get that for any $\kappa$ $$\sup_{n=0,\dots,N}\|U^n_{\kappa} -U^n\|_\infty\le \frac{\Gamma(2-\alpha)}{\alpha(1-\alpha)}T^\alpha    E_{\kappa}$$ and therefore the continuity of the map $\Psi_1$.\
Given $U=\Psi_1(M)$, we first show that, for ${{\Delta t}}$ sufficiently small, $\tilde M=\Psi_2(U)$ is well defined. First observe that, by Prop. \[prop:lip\_U\], $\|D_h U\|_\infty$ is bounded uniformly in $M$. Therefore, since $g$ is $C^1$, there exists a constant $C$, independent of $M$, such that for $k=1,2,3,4,$ $$\label{eq:Fix1}
\left|\frac{\partial g}{\partial q_k} (x_{i,j},[D_hU^n]_{i,j})\right|\le C, \quad 0\le i,j\le N_h,\,0\le n\le N.$$ Hence, arguing as in Prop. \[prop:ex\_FP\], we obtain that, for ${{\Delta t}}$ sufficiently small, there exists a solution $\tilde M$ of . Moreover $\tilde M\in {\mathcal{K}}_h^{N}$. The continuity of the map $\Psi_2$ follows by the regularity of $g$ and the linearity of the problem.\
The boundedness and continuity of the map $\Phi_1$ and the continuity of the map $\Psi_2$ implies that $\Phi=\Psi_2\circ\Psi_1$ is a continuous map on the compact, convex set ${\mathcal{K}}_h^N$. Moreover by Prop. \[prop:ex\_FP\], $\Phi$ maps ${\mathcal{K}}_h^N$ into itself. Hence, by Brouwer’s Theorem, $\Phi$ admits a fixed point $M$, which clearly a solution of system .

Convergence {#sec:convergence}
===========

In this section we prove a convergence result for the scheme . We relies on the analysis performed in [@acc Section 3] for the Mean Field Games system and in [@lx Theorem 3.2] for computation related to fractional derivatives.\
A preliminary result for the analysis is a fundamental identity which relies on the duality structure of the discrete system.

Let $(U,M)$ be a solution of and $({\tilde U},{\tilde M})$ a solution of the perturbed system $$\label{eq:conv_scheme_pert}
\left\{
\begin{array}{ll}
\bar D_{\Delta t}^{\alpha} {\tilde U}_{i,j}^{n} - (\Delta_h {\tilde U}^{n})_{i,j} + g(x_{i,j},[D_h{\tilde U}^{n}]_{i,j}) = f_h[{\tilde M}^{n+1}]_{i,j}+a^n_{i,j}, \\[6pt]
D_{\Delta t}^\alpha {\tilde M}_{i,j}^{n+1} - (\Delta_h {\tilde M}^{n+1})_{i,j} - {\mathcal{B}}_{i,j}({\tilde U}^{n},{\tilde M}^{n+1}) = b^{n+1}_{i,j}, \\[6pt]
{\tilde U}^{N}_{i,j} = u_T(x_{i,j}), \,   {\tilde M}_{i,j}^{0}=  \frac 1 {h^2} \int_{|x-x_{i,j}|_{\infty}\le h/2} m_0(x) dx,
\end{array}
\right.$$ for $n=0,\dots,N-1$ and $0\le i,j\le N_h$. Then $$\label{eq:conv_fund_id}
\begin{aligned}
&  \sum_{n=0}^{N-1}   \bigg(f_h[M^{n+1}] - f_h[\tilde{M}^{n+1}],M^{n+1} - \tilde{M}^{n+1}\bigg)_2 \\
& +   \sum_{n=0}^{N-1}  \sum_{i,j} M_{i,j}^{n+1}R^n_{i,j}(U,{\tilde U})+ \sum_{n=0}^{N-1}  \sum_{i,j} \tilde{M}_{i,j}^{n+1} R^n_{i,j}({\tilde U},U)\\
&=\sum_{n=0}^{N-1}(a^n,M^{n+1}-{\tilde M}^{n+1})_2+\sum_{n=0}^{N-1}(b^{n+1},U^n-{\tilde U}^n)_2
\end{aligned}$$ where $$R^n_{i,j}(V,W)=g(x_{i,j},[D_hW^{n}]_{i,j}) - g(x_{i,j},[D_hV^{n}]_{i,j})   -[D_h(W^{n}-V^{n})]_{i,j} \cdot  g_q  (x_{i,j},[D_hV^{n}])$$ for $V=(V^n)_{0\le n\le N}$, $W=(W^n)_{0\le n\le N}$.

The identity can be obtained by arguing as in the proof of Theorem \[thm:uniqueness\], see in particular , and taking into account the additional terms $a^n_{i,j}$, $b^n_{i,j}$.

For the convergence analysis, we consider the reference case $$\label{eq:conv_ham}
H(x,p)=|p|^\beta+c(x)$$ where $\beta\ge 2$ and $c:{\mathbb{T}}^2\to{\mathbb{R}}$ is a continuous function, but the analysis can be extended in a more general setting. For the approximation of the Hamiltonian, we consider the numerical Hamiltonian $$\label{eq:conv_num_ham}
g(x,q)=((q_1^-)^2+(q_2^+)^2+(q_3^-)^2+(q_4^+)^2)^{\beta/2}+c(x)$$ where $q^\pm$ denote the positive and the negative part of $q\in{\mathbb{R}}$.

Assume that the continuous problem admits a unique classical solution $(u,m)\in C^{2}({\mathbb{T}}^2\times [0,T])\times C^{2}({\mathbb{T}}^2\times [0,T])$ and that $m_0\ge \delta>0$. For ${\sigma}=(h,{{\Delta t}})$, denote by $(u_{\sigma},m_{\sigma})$ the linear function obtained by interpolating the solution $(U,M)$ of on the nodes of the space-time grids. Then, as $|{\sigma}|\to 0$, $(u_{\sigma},m_{\sigma})$ converges to $(u,m)$ in $L^\beta(0,T: W^{1,\beta}({\mathbb{T}}^2))\times  (C(0,T;L^2({\mathbb{T}}^2)) \cap L^2(0,T;H^1({\mathbb{T}}^2))$.

The idea of the proof is to show that the solution of the continuous problem satisfies the perturbed discrete problem with $(a^n,b^n)$ vanishing for the discretion steps ${{\Delta t}},h$ going to $0$.\
Given a solution $(u,m)$ of , we define two grid functions ${\tilde U}_{\sigma}$,${\tilde M}_{\sigma}$ on $({\mathbb{R}}^{N_h^2})^N$ by $$({\tilde U}_{\sigma})^n_{i,j}=u(x_{i,j},n{{\Delta t}}),\quad ({\tilde M}_{\sigma})^n_{i,j}=\frac 1 {h^2} \int_{|x-x_{i,j}|\le h/2} m(x,n{{\Delta t}}) dx .$$ Note that ${\tilde M}_{\sigma}\in {\mathcal{K}}_h^N$, see . Moreover, since $(u,m)$ is a smooth function, $({\tilde U}_{\sigma},{\tilde M}_{\sigma})$ is a solution of with the perturbation $(a^n_{\sigma},b^n_{\sigma})$ given by the consistency error of the scheme. Hence $$\label{eq:conv_consistency}
\max_{0\le n\le N}\{\|a^n_{\sigma}\|_\infty+\|b^n_{\sigma}\|_\infty\}\to 0 \quad \text{for ${\sigma}\to 0$}.$$ The first step of the proof is the convergence result $$\label{eq:conv_seminorm}
\lim_{|{\sigma}|\to 0}h^2{{\Delta t}}\sum_{n=0}^N\sum_{i,j}|[D_hU^n_{\sigma}]_{i,j}-[D_h({\tilde U}_{\sigma})^n]_{i,j}|^\beta=0.$$ The identity is obtained by means of the fundamental identity which does not involve the discrete fractional derivatives. Hence, for its proof, it is possible to repeat the same proof of Step 1 in [@acc Section 3]. Note that here it is employed that ${\tilde M}_{\sigma}>0$, which follows by the assumption $m_0\ge\delta>0$.\
The second step of the proof is the convergence of the solution of the Fokker-Planck equation. We denote with $M$, ${\tilde M}$ the grid functions corresponding to $m_{\sigma}$ and $m$, omitting the index ${\sigma}$ for simplicity. We set $E=M-{\tilde M}$ and $$\begin{aligned}
&\|[D_hE^{n+1}]\|^2_2=h^2\sum_{i,j} [D_hE^{n+1}]_{i,j}\cdot[D_hE^{n+1}]_{i,j},\\
&\|(E^{n+1},[D_hE^{n+1}])\|_2^2=\|E^{n+1}\|^2_2+\rho_\alpha\|[D_hE^{n+1}]\|^2_2.\end{aligned}$$ By subtracting the equations satisfied by $M$, ${\tilde M}$, we get $$E^{n+1}_{i,j}-\rho_\alpha (\Delta_h E^{n+1})_{i,j} - \rho_\alpha({\mathcal{B}}_{i,j}(U^{n},M^{n+1})-{\mathcal{B}}_{i,j}({\tilde U}^{n},{\tilde M}^{n+1})) =\sum_{k=0}^n c^n_k E^k+{\rho_\alpha}b^{n+1}.$$ Multiplying the previous equation by $E^{n+1}_{i,j}$, summing over $i$ and $j$ and using the identities and , we get $$\label{eq:conv_es1}
\begin{split}
&\|(E^{n+1},[D_hE^{n+1}])\|_2^2-{\rho_\alpha}h^2 \sum_{i,j} E^{n+1}_{i,j} [D_hE^{n+1}]_{i,j}\cdot g_q(x_{i,j}, [D U^n]_{i,j})\\
&=h^2\sum_{k=0}^n c^n_k (E^k,E^{n+1})+{\rho_\alpha}h^2(b^{n+1}, E^{n+1})\\
&+{\rho_\alpha}h^2  \sum_{i,j}{\tilde M}_{i,j}^{n+1}[DE^{n+1}]_{i,j}\cdot(g_q(x_{i,j}, [D_hU^n]_{i,j})-g_q(x_{i,j}, [D_h{\tilde U}^n]_{i,j})).
\end{split}$$ Taking into account , we estimate for ${{\Delta t}}$ sufficiently small $$\left|{\rho_\alpha}h^2 \sum_{i,j} E^{n+1}_{i,j} [D_hE^{n+1}]_{i,j}\cdot g_q(x_{i,j}, [D U^n]_{i,j})\right|\le 
\frac 1 2\|(E^{n+1},[D_hE^{n+1}])\|_2^2.$$ Using again and since $(u,m)$ is a classical solution to , we have for all $\eta>0$ (see[@acc Equation (4.5)]) $$\begin{aligned}
&\left|{\tilde M}_{i,j}^{n+1}[DE^{n+1}]_{i,j}\cdot(g_q(x_{i,j}, [D_hU^n]_{i,j})-g_q(x_{i,j}, [D_h{\tilde U}^n]_{i,j}))\right|\\
&\le  C \Big(\frac{1}{\eta} |[D_hE^{n+1}]_{i,j}|^2+\eta |[D_hU^n]_{i,j}-[D_h{\tilde U}^n]_{i,j}|^\beta\Big)\end{aligned}$$ Plugging the previous estimates in and dividing by $\|(E^{n+1},[D_hE^{n+1}])\|_2^2$, we get $$\label{eq:conv_es2}
\|(E^{n+1},[D_hE^{n+1}])\|_2\le C\left(\sum_{k=0}^n c^n_k \|E^k\|_2+{\rho_\alpha}\varepsilon^{n+1}\right)$$ where $$\varepsilon^{n+1}=\|b^{n+1}\|_\infty+  h^2  \sum_{i,j}|[D_hU^n]_{i,j}-[D_h{\tilde U}^n]_{i,j}|^\beta.$$ and the constant $C$ is independent of $h,{{\Delta t}}$.\
Following [@lx Theorem 3.2], we rewrite as $$\label{eq:conv_es3}
\|(E^{n+1},[D_hE^{n+1}])\|_2\le C\Big((1-B_1)\|E^n\|_2+\sum_{k=1}^{n-1}(B_k-B_{k+1})\|E^{n-k}\|_2
+B_n\|E^0\|_2+{\rho_\alpha}\varepsilon^{n+1}\Big),$$ where $B_k=(k+1)^{1-\alpha}-k^{1-\alpha}$, $k=0,\dots,n$. Recall that, since $M^0={\tilde M}^0$, then $E^0\equiv 0$. We claim that for $n=1,\dots,N$, $$\label{eq:conv_es4}
\|(E^{n},[D_hE^{n}])\|_2\le C{\rho_\alpha}B_{n-1}\max_{0\le k\le N}\varepsilon^k.$$ For $n=1$, by we have $$\|(E^{1},[D_hE^{1}])\|_2\le C{\rho_\alpha}\varepsilon^1\le C {\rho_\alpha}B_0^{-1}\max_{0\le k\le N}\varepsilon^k,$$ since $B_0=1$. Assuming that is true for $n$, by and the fact that $B_k^{-1}/B_{k+1}^{-1}<1$, we obtain $$\begin{aligned}
&\|(E^{n+1},[D_hE^{n+1}])\|_2\le \left[ (1-B_1)B_k^{-1}+\sum_{k=1}^{n-1}(B_k-B_{k+1})B^{-1}_{n-k-1}\right]C{\rho_\alpha}\max_{0\le k\le N}\varepsilon^k\\
&+C{\rho_\alpha}\max_{0\le k\le N}\varepsilon^k
\le \left[ (1-B_1) +\sum_{k=1}^{n-1}(B_k-B_{k+1})+B_n\right]C  {\rho_\alpha}B_n^{-1}\max_{0\le k\le N}\varepsilon^k.\end{aligned}$$ Using $(1-B_1) +\sum_{k=1}^{n-1}(B_k-B_{k+1})+B_n=1$ (see [@lx Equation (3.7)]), the claim follows.\
Since $n^{-\alpha}B_{n-1}^{-1}\le 1/(1-\alpha)$ for $n=1,\dots,N$ (see [@lx Equation (3.21)]) and ${\rho_\alpha}=\Gamma(2-\alpha)\Delta t^{\alpha}$, by we have $$\label{eq:conv_es5}
\|(E^{n},[D_hE^{n}])\|_2\le C\max_{0\le k\le N}\varepsilon^k.$$ Recalling that $E^n$ and $\varepsilon^n$ depends on ${\sigma}$ and that, by -, $\lim_{|{\sigma}|\to 0}\max_{0\le k\le N}\varepsilon^k_{\sigma}=0$, we finally get $$\lim_{|{\sigma}|\to 0}\,\max_{0\le n\le N}\{\|E_{\sigma}^{n}\|_2+\rho_\alpha\|[D_hE_{\sigma}^{n}]\|_2\}=0,$$ which gives the convergence of $m_{\sigma}$ to $m$ in $C(0,T;L^2({\mathbb{T}}^2))$. To get the convergence in $L^2(0,T;H^1({\mathbb{T}}^2))$, it is sufficient to observe that, since $m$ is smooth, the consistency error of the approximation of the fractional derivative is of order ${{\Delta t}}^{2-\alpha}$ (see [@lx Equation (3.3)]). Hence, by , it follows $\lim_{|{\sigma}|\to 0} {{\Delta t}}\sum_{n=1}^N\|[D_hE_s^{n}]\|_2=0$ and therefore the convergence of the gradient of the error.\
In the final step of the proof, we prove the convergence of $u_{\sigma}$ to $u$. By the convergence of $m_{\sigma}$ to $m$ and assumption (F4), we have that ${\tilde U}$ is a solution to a discrete Hamilton-Jacobi-Bellman equation with $M$ in place of ${\tilde M}$, i.e. $$\label{eq:conv_HJ_pert}
\left\{
\begin{array}{ll}
\bar D_{\Delta t}^{\alpha} {\tilde U}_{i,j}^{n} - (\Delta_h {\tilde U}^{n})_{i,j} + g(x_{i,j},[D_h{\tilde U}^{n}]_{i,j}) = f_h[M^{n+1}]_{i,j}+a^n_{i,j},\\[6pt] 
{\tilde U}^{N_T}_{i,j} = u_T(x_{i,j}),
\end{array}
\right.$$ where the perturbation $a^n$ still satisfies . Arguing as in Prop. \[prop:lip\_U\] and with the same notation, we have $U^n= \Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k U^k\right)$, ${\tilde U}^n= \Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k {\tilde U}^k+a^n\right)$. Hence $$\begin{aligned}
&\|{\tilde U}^n-U^n\|_\infty= \left\|\Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k {\tilde U}^k+\rho_\alpha a^n\right)
-\Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k  U^k\right)\right\|_\infty\\
&\le \sum_{k=n+1}^{N}{\overline c}_{n}^k \left\| {\tilde U}^k - U^k\right\|_\infty+{\rho_\alpha}\|a^n\|_\infty.\end{aligned}$$ The previous inequality is equivalent to $$\bar D_{\Delta t}^{\alpha}(\|{\tilde U}^n-U^n\|_\infty)\le \|a^n\|_\infty.$$ We conclude, by the backward version of discrete Gronwall inequality , the uniform convergence of $u_{\sigma}$ to $u$. The convergence of the gradients follows by .

In this paper we have focused on the numerical approximation of the classical solution of the PDE system and we have assumed the existence of a unique classical solution. We refer to [@q] for the existence and uniqueness of a weak solution to the system, where the solution of the fractional Fokker-Planck equation is defined in a time-fractional Bochner space and the equation is considered in the sense of distributions. It would be very interesting to define some notion of renormalized solution to fractional PDEs and consider weak solution to the system, as well as approximation schemes. In the standard case, convergence of a finite differences scheme to the weak solution of a Mean Field Games system was shown in [@ap].

Numerical simulations {#sec:numerics}
=====================

In this section, we illustrate some numerical simulations in one-dimensional case. In all examples, we consider the domain ${\mathbb{T}}= [0,1]$ and the final time $T=2$. We test $$\label{eq:test}
\begin{cases}
\partial_{[t,T)}^\alpha u -\sigma \Delta u + \frac{u_x^2}{2} = F(x,t,m) \\[4pt]
\partial_{(0,t]}^\alpha m - \sigma \Delta m - (mu_x)_x = 0 \\[4pt]
\end{cases}$$ with terminal-initial data $u(x,T)=0$ and $m(x,0)=\nu(x)/\int_{\mathbb{T}}\nu dx$ where $\nu(x) = e^{-(x-0.5)^2/(0.1)^2}$. Here, the running cost function has a form $$F(x,t,m) = 5\left(x-\frac{1-\sin (2 \pi t)}{2}\right)^2 + \lambda m$$ where $\lambda$ is a parameter.

**Test 1.** For $(x,t)\in{\mathbb{T}}\times [0,2]$, we take $\sigma = 0$ and $\lambda = 0$, i.e. there is no diffusion and penalization term in . Even if this case is not covered by the results of the paper, we consider it to highlight the effects of the fractional derivative. For the running cost in the MFG system, the agents try to concentrate at the point $(1-\sin (2 \pi t))/2$. We apply via full discretization of . We take $N_h=50$ and $N=2000$ so that $h=\frac{1}{50}$ and $\Delta t = \frac{T}{2000}$. Since in this case the Hamilton-Jacobi-Bellman equation in the system does not depend on $m$, the system is decoupled. In the numerical computation, we start from the terminal data $U_i^{N_T} = u_T(x_i)$ and go backward to iteratively calculate $U^n$ on the whole grid for all $0 \le n < N_T$. Having the numerical solution for the Hamilton-Jacobi-Bellman equation, we now can compute $M^n$ on the grid for all $0<n \le N$ starting from the initial time. In Figure 1, we provide density evolution for three different values of $\alpha$. From the level sets in Figure 2, we can see that the fractional time derivative has a kind of diffusion effect. Also, from Figure 2, we notice that when $\alpha$ gets smaller, the shape of the sinusoidal curve becomes more and more asymmetric and smoother. This means that for the smaller values of $\alpha$, in the mean-field game, the agents try to spread away from the sinusoidal curve.

  ---------------------------------------------------- ------------------------------------------------------ -----------------------------------------------------
   ![image](figures/test1_m_alpha1.png){width="50mm"}   ![image](figures/test1_m_alpha085.png){width="50mm"}   ![image](figures/test1_m_alpha07.png){width="50mm"}
                          (A)                                                   (B)                                                    (C)
  ---------------------------------------------------- ------------------------------------------------------ -----------------------------------------------------

\

  -------------------------------------------------------- ---------------------------------------------------------- ---------------------------------------------------------
   ![image](figures/test1_m_alpha1_ls1.png){width="50mm"}   ![image](figures/test1_m_alpha1_ls085.png){width="50mm"}   ![image](figures/test1_m_alpha1_ls07.png){width="50mm"}
                            (A)                                                       (B)                                                        (C)
  -------------------------------------------------------- ---------------------------------------------------------- ---------------------------------------------------------

\

**Test 2.** Now, we consider with $\sigma = 0.1$ and $\lambda=0$ to see how the diffusion term affects the density evolution. For this case, we provide only level sets of the density in Figure 3. During the whole time interval, a diffusive effect is observed. Moreover, its effect is much stronger than the fractional time derivative effect which is displayed in Figure 2. For example, we can compare Figure 2 (B), where $\alpha=0.85, \sigma =0$ and Figure 3 (A), where $\alpha=1, \sigma =0.1$.

  ----------------------------------------------------------- ------------------------------------------------------------- ------------------------------------------------------------
   ![image](figures/test2_levelset_alpha1.png){width="50mm"}   ![image](figures/test2_levelset_alpha085.png){width="50mm"}   ![image](figures/test2_levelset_alpha07.png){width="50mm"}
                              (A)                                                          (B)                                                          (C)
  ----------------------------------------------------------- ------------------------------------------------------------- ------------------------------------------------------------

\

**Test 3.** We test the case where $\sigma=0$ and $\lambda=1$. Since the coupling cost is increasing in $m$, it penalizes concentration of the agents. Indeed, we see from Figure 4 that concentration of agents is not as high as in Test 1.

  ----------------------------------------------------------- ------------------------------------------------------------- ------------------------------------------------------------
   ![image](figures/test3_levelset_alpha1.png){width="50mm"}   ![image](figures/test3_levelset_alpha085.png){width="50mm"}   ![image](figures/test3_levelset_alpha07.png){width="50mm"}
                              (A)                                                          (B)                                                          (C)
  ----------------------------------------------------------- ------------------------------------------------------------- ------------------------------------------------------------

\

Discussion
==========

In this paper, we considered the numerical approximation to a time-fractional forward backward PDE system. By using a classical scheme from Mean Field Games theory and L1 scheme for fractional derivatives the duality structure of the PDE system is preserved on the discrete level. The duality with regard to the fractional derivatives is further clarified in the following discrete fractional integration by parts formula in the appendix. In future work we may consider the error estimates of approximation scheme to the viscosity solution of the fractional Hamilton-Jacobi-Bellman equation as in [@glm]. Another interesting direction will be to consider solving the fractional forward backward system by machine learning methods. Machine learning methods for solving mean field games have been considered by Lauriere and Carmona in [@ml]. In [@Pang], Pang, Li and Karniadakis used machine learning methods to solve forward and inverse problems of time-space fractional PDEs. A combination of techniques might produce fruitful results for studying fractional PDE systems.

Appendix
========

In this proof, we use (ii) in Lemma \[lem:c\_properties\]. We have $$\begin{aligned}
{\rho_\alpha}& \sum_{n=0}^{N-1} \left( \bar D_{\Delta t}^{\alpha} U^{n}, M^{n+1} \right)_2 + \sum_{n=0}^{N-1} \bar c_n^N \left(U^N,M^{n+1}\right)_2  \\
& = \sum_{n=0}^{N-1} \left(U^n - \sum_{k=n+1}^{N} \bar c_n^k U^k,M^{n+1}\right)_2 + \sum_{n=0}^{N-1} \bar c_n^N \left(U^N,M^{n+1}\right)_2 \\
& = \sum_{n=0}^{N-1} (U^n,M^{n+1})_2 - \sum_{n=0}^{N-2} \sum_{k=n+1}^{N-1} \left( \bar c_n^k U^k,M^{n+1}\right)_2 \quad \text{(change the order of sums)}\\
& =  \sum_{n=0}^{N-1} (U^n,M^{n+1})_2 - \sum_{k=1}^{N-1} \sum_{n=0}^{k-1} \left( \bar c_n^k U^k,M^{n+1}\right)_2 \quad\quad \text{(interchange $m$ and $n$)}\\
& = \sum_{n=0}^{N-1} (U^n,M^{n+1})_2 - \sum_{n=1}^{N-1} \sum_{k=0}^{n-1} \left( \bar c_k^n U^n,M^{k+1}\right)_2\\
& = (U^0,M^1)_2 + \sum_{n=1}^{N-1} (U^n, M^{n+1} - \sum_{k=0}^{n-1} \bar c_{k}^{n} M^{k+1})_2 \\
& = (U^0,M^1)_2 + \sum_{n=1}^{N-1} (U^n, M^{n+1} - \sum_{k=0}^{n-1} c_{k+1}^{n+1} M^{k+1})_2 \quad\quad  \text{(since $\bar c_{k}^{n} = c_{k+1}^{n+1}$)}\\
& = (U^0,M^1)_2 + \sum_{n=1}^{N-1} (U^n, M^{n+1} - \sum_{k=1}^{n} c_{k}^{n+1} M^{k})_2 \\
& = (U^0,M^1 - c_0^1 M^0)_2 + \sum_{n=1}^{N-1} (U^n, M^{n+1} - \sum_{k=0}^{n} c_{k}^{n+1} M^{k})+\sum_{n=0}^{N-1} c_0^{n+1}(U^n,M^0)_2 \\
& = \sum_{n=0}^{N-1} (U^n, M^{n+1} - \sum_{k=0}^{n} c_{k}^{n+1} M^{k})_2 + \sum_{n=0}^{N-1} c_0^{n+1}(U^n,M^0)_2 \\
& = {\rho_\alpha}\sum_{n=0}^{N-1} (D_{\Delta t}^{\alpha} M^{n+1},U^{n})_2 + \sum_{n=0}^{N-1} c_0^{n+1}(U^n,M^0)_2.
\end{aligned}$$ Dividing by ${\rho_\alpha}$, we obtain the result.

[99]{}

Achdou, Y.; Capuzzo Dolcetta, I. Mean field games: numerical methods. SIAM J. Numer. Anal. 48 (2010), no. 3, 1136–1162.

Achdou, Y.; Camilli, F.; Capuzzo Dolcetta, I. Mean field games: convergence of a finite difference method. SIAM J. Numer. Anal. 51 (2013), no. 5, 2585-2612.

Achdou, Y.; Porretta, A. Convergence of a finite difference scheme to weak solutions of the system of partial differential equations arising in mean field games. SIAM J. Numer. Anal. 54 (2016), 161-186.

Almulla, N.; Ferreira, R.; Gomes, D. Two numerical approaches to stationary mean-field games. Dyn. Games Appl. 7 (2017), no. 4, 657-682.

Annunziato, M.; Borzì, A.; Magdziarz, M.; Weron, A. A fractional Fokker-Planck control framework for subdiffusion processes. Optimal Control Appl. Methods 37 (2016), no. 2, 290-304.

Bouchaud, J.-P.; Georges, A. Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications. Phys. Rep. 195 (1990), no. 4-5, 127-293.

Briceño-Arias, L. M.; Kalise, D.; Silva, F. J. Proximal methods for stationary mean field games with local couplings. SIAM J. Control Optim. 56 (2018), no. 2, 801-836.

Camilli, F.; Goffi, A. Existence and regularity results for viscous Hamilton-Jacobi equations with Caputo time-fractional derivative, NoDEA Nonlinear Differential Equations Appl. 27 (2020), no. 2, Paper No. 22.

Carlini, E.; Silva, F. A semi-Lagrangian scheme for a degenerate second order mean field game system. Discrete Contin. Dyn. Syst. 35 (2015), no. 9, 4269-4292.

Du, Q.; Yang, J.; Zhou, Z. Time-fractional Allen-Cahn equations: Analysis and numerical methods. , 1906.06584, 2019.

Du, Q. An invitation to nonlocal modeling, analysis and computation. , 2018, volume 3, 3523–3552.

Giga, Y; Liu, Q.; Mitake, H. On a discrete scheme for time fractional fully nonlinear evolution equations, arXiv:1904.06105.

Gorenflo, R.; Luchko, R.; Yamamoto, M. Time-fractional diffusion equation in the fractional Sobolev spaces. Fract. Calc. Appl. Anal. 18 (2015), no. 3, 799–820.

Jin, B.; Lazarov, R.; Pasciak, J; Rundell, W. Variational formulation of problems involving fractional order differential operators. Math. Comp. 84 (2015), no. 296, 2665–2700.

Jin, B.; Lazarov, R.; Zhou, Z. An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal. 36 (2016), no. 1, 197-221.

Jin, B.; Lazarov, R.; Zhou, Z. Numerical methods for time-fractional evolution equations with nonsmooth data: A concise overview. Comput. Methods Appl. Mech. Engrg. 346 (2019), 332–358.

Kolokoltsov, V; Veretennikova, M. A fractional Hamilton-Jacobi Bellman equation for scaled limits of controlled continuous time random walks. Commun. Appl. Ind. Math. 6 (2014), no. 1, e-484, 18 pp.

Lasry, J.-M.; Lions, P.-L. Mean Field Games. Jpn. J. Math., 2 (2007), 229-260.

Lauriere, M; Carmona, R. Convergence Analysis of Machine Learning Algorithms for the Numerical Solution of Mean Field Control and Games: II – The Finite Horizon Case. arXiv:1908.01613v1, 2019.

Li, D.; Wang, J.; Zhang, J. Unconditionally convergent L1-Galerkin FEMs for nonlinear time-fractional Schrödinger equations. SIAM J. Sci. Comput. 39 (2017), no. 6, A3067-A3088.

Lin, Y.; Xu, C. Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput. Phys. 225 (2007), no. 2, 1533–1552.

Metzler, R.; Klafter, J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A 37 (2004), no. 31, R161-R208.

Meerschaert, M.; Sikorskii, A. [*Stochastic Models for Fractional Calculus*]{}, De Gruyter Studies in Mathematics, 43. De Gruyter, Berlin, 2019.

Namba, T. On existence and uniqueness of viscosity solutions for second order fully nonlinear pdes with caputo time fractional derivatives. , 25 (2018), no. 3, Paper No. 23,

Pang, G.; Lu, L.; Karniadakis, G. fpinns: Fractional physics-informed neural networks. , 41 (2019), A2603–A2626.

Shen, J.; Sheng, C.-T. An efficient space–time method for time fractional diffusion equation. , 81 (2019), 1088–1110.

Tang, Q. On an optimal control problem of time-fractional advection-diffusion equation. Discrete & Continuous Dynamical Systems-B, 25 (2020), 761-779.

Tang, T.; Yu, H; Zhou., T. On energy dissipation theory and numerical stability for time-fractional phase-field equations. , 41 (2019), A3757–A3778.

Topp, E.; Yangari, M.. Existence and uniqueness for parabolic problems with Caputo time derivative. , 262 (2017), 6018 – 6046.

Zaslavsky, G. Chaos, fractional kinetics, and anomalous transport. , 371 (2002), 461–580.
