  We propose a numerical method for the Vlasov-Poisson-Fokker-Planck model
written as an hyperbolic system thanks to a spectral decomposition in the basis
of Hermite functions with respect to the velocity variable and a structure
preserving finite volume scheme for the space variable. On the one hand, we
show that this scheme naturally preserves both stationary solutions and
linearized free-energy estimate. On the other hand, we adapt previous arguments
based on hypocoercivity methods to get quantitative estimates ensuring the
exponential relaxation to equilibrium of the discrete solution for the
linearized Vlasov-Poisson-Fokker-Planck system, uniformly with respect to both
scaling and discretization parameters. Finally, we perform substantial
numerical simulations for the nonlinear system to illustrate the efficiency of
this approach for a large variety of collisional regimes (plasma echos for
weakly collisional regimes and trend to equilibrium for collisional plasmas)
and to highlight its robustness (unconditional stability, asymptotic preserving
properties).
