In this work, we present efficient and high-order numerical schemes for sub-diffusion equation involving tempered fractional derivative. In particular, we establish the equivalence between the sub-diffusion equation with tempered fractional derivative in a d dimensional space and integer-order extended parametric sub-diffusion equation in a \((d+1)\) dimensional space, following our previous work [J. Sci. Comput.,105(2025)82]. We apply the BDF-k scheme for time discretization and conduct a rigorous stability analysis for the semi-discrete schemes. Motivated by the high-regularity with respect to the extended \(\theta \) direction established in [J. Sci. Comput.,105(2025)82], we employ a Jacobi spectral method for the extended \(\theta \) direction and Fourier spectral method for the spatial direction, respectively. Additionally, we provide an error estimate for the full discretization scheme, showing a convergence order of \(\mathcal {O}\left( \Delta t^{k} + M^{-m} + N^{-q}\right) , \, 1\le k \le 5, \, 0<m\le M\) , where \(\Delta t\) , M, and N represent the time-step size, the number of nodes for the Jacobi spectral method and Fourier spectral method, respectively, q is determined by the regularity of the solution in the spatial direction. To accelerate the computation, we apply a generalized characteristic decomposition to the fully discrete scheme, which yields M independent discrete Poisson systems. The computational cost and storage requirements of our proposed algorithm are essentially the same as those for discrete Poisson systems, namely, \(\mathcal {O}\left( N_TN^2\log (N^2)\right) \) in cost and \(\mathcal {O}\left( N^2\right) \) in storage, where \(N_T\) is the number of time-step. We present several numerical examples to demonstrate the effectiveness of the proposed method.
