  In this paper, we consider a finite difference grid-based semi-Lagrangian
approach in solving the Vlasov-Poisson (VP) system. Many of existing methods
are based on dimensional splitting, which decouples the problem into solving
linear advection problems, see {\em Cheng and Knorr, Journal of Computational
Physics, 22(1976)}. However, such splitting is subject to the splitting error.
If we consider multi-dimensional problems without splitting, difficulty arises
in tracing characteristics with high order accuracy. Specifically, the
evolution of characteristics is subject to the electric field which is
determined globally from the distribution of particle densities via the
Poisson's equation. In this paper, we propose a novel strategy of tracing
characteristics high order in time via a two-stage multi-derivative
prediction-correction approach and by using moment equations of the VP system.
With the foot of characteristics being accurately located, we proposed to use
weighted essentially non-oscillatory (WENO) interpolation to recover function
values between grid points, therefore to update solutions at the next time
level. The proposed algorithm does not have time step restriction as Eulerian
approach and enjoys high order spatial and temporal accuracy. However, such
finite difference algorithm does not enjoy mass conservation; we discuss one
possible way of resolving such issue and its potential challenge in numerical
stability. The performance of the proposed schemes are numerically demonstrated
via classical test problems such as Landau damping and two stream
instabilities.
