  We propose and analyze a mixed finite element method for the spatial
approximation of a time-fractional Fokker--Planck equation in a convex
polyhedral domain, where the given driving force is a function of space. Taking
into account the limited smoothing properties of the model, and considering an
appropriate splitting of the errors, we employed a sequence of clever energy
arguments to show optimal convergence rates with respect to both approximation
properties and regularity results. In particular, error bounds for both primary
and secondary variables are derived in $L^2$-norm for cases with smooth and
nonsmooth initial data. We further investigate a fully implicit time-stepping
scheme based on a convolution quadrature in time generated by the backward
Euler method. Our main result provides pointwise-in-time optimal $L^2$-error
estimates for the primary variable. Numerical examples are then presented to
illustrate the theoretical contributions.
