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. 