  We design a variational asymptotic preserving scheme for the
Vlasov-Poisson-Fokker-Planck system with the high field scaling, which
describes the Brownian motion of a large system of particles in a surrounding
bath. Our scheme builds on an implicit-explicit framework, wherein the stiff
terms coming from the collision and field effects are solved implicitly while
the convection terms are solved explicitly. To treat the implicit part, we
propose a variational approach by viewing it as a Wasserstein gradient flow of
the relative entropy, and solve it via a proximal quasi-Newton method. In so
doing we get positivity and asymptotic preservation for free. The method is
also massively parallelizable and thus suitable for high dimensional problems.
We further show that the convergence of our implicit solver is uniform across
different scales. A suite of numerical examples are presented at the end to
validate the performance of the proposed scheme.
