We provide and analyze the high order algorithms for the model describing the functional distributions of particles performing anomalous motion with power-law jump length and tempered power-law waiting time. The model is derived in Wu et al. (Phys Rev E 93:032151, 2016), being called the time-tempered fractional Feynman–Kac equation named after Richard Feynman and Mark Kac who first considered the model describing the functional distribution of normal motion. The key step of designing the algorithms is to discretize the time tempered fractional substantial derivative, being defined as SDtγ,λ~G(x,p,t)=Dtγ,λ~G(x,p,t)-λγG(x,p,t)\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\begin{aligned} {^S\!}D_t^{\gamma ,\widetilde{\lambda }} G(x,p,t)\!=\!D_t^{\gamma ,\widetilde{\lambda }} G(x,p,t)\!-\!\lambda ^\gamma G(x,p,t) \end{aligned}$$\end{document}with λ~=λ+pU(x),p=ρ+Jη,J=-1\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\widetilde{\lambda }=\lambda + pU(x),\, p=\rho +J\eta ,\, J=\sqrt{-1}$$\end{document}, where Dtγ,λ~G(x,p,t)=1Γ(1-γ)∂∂t+λ~∫0tt-z-γe-λ~·(t-z)G(x,p,z)dz,\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\begin{aligned} D_t^{\gamma ,\widetilde{\lambda }} G(x,p,t) =\frac{1}{\varGamma (1-\gamma )} \left[ \frac{\partial }{\partial t}+\widetilde{\lambda } \right] \int _{0}^t{\left( t-z\right) ^{-\gamma }}e^{-\widetilde{\lambda }\cdot (t-z)}{G(x,p,z)}dz, \end{aligned}$$\end{document}and λ≥0\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\lambda \ge 0$$\end{document}, 0<γ<1\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$0<\gamma <1$$\end{document}, ρ>0\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\rho >0$$\end{document}, and η\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\eta $$\end{document} is a real number. The designed schemes are unconditionally stable and have the global truncation error O(τ2+h2)\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\mathscr {O}(\tau ^2+h^2)$$\end{document}, being theoretically proved and numerically verified in complex space. Moreover, some simulations for the distributions of the first passage time are performed, and the second order convergence is also obtained for solving the ‘physical’ equation (without artificial source term).
