  In this paper we propose and analyze a finite difference numerical scheme for
the Poisson-Nernst-Planck equation (PNP) system. To understand the energy
structure of the PNP model, we make use of the Energetic Variational Approach
(EnVarA), so that the PNP system could be reformulated as a non-constant
mobility $H^{-1}$ gradient flow, with singular logarithmic energy potentials
involved. To ensure the unique solvability and energy stability, the mobility
function is explicitly treated, while both the logarithmic and the electric
potential diffusion terms are treated implicitly, due to the convex nature of
these two energy functional parts. The positivity-preserving property for both
concentrations, $n$ and $p$, is established at a theoretical level. This is
based on the subtle fact that the singular nature of the logarithmic term
around the value of $0$ prevents the numerical solution reaching the singular
value, so that the numerical scheme is always well-defined. In addition, an
optimal rate convergence analysis is provided in this work, in which many
highly non-standard estimates have to be involved, due to the nonlinear
parabolic coefficients. The higher order asymptotic expansion (up to third
order temporal accuracy and fourth order spatial accuracy), the rough error
estimate (to establish the $\ell^\infty$ bound for $n$ and $p$), and the
refined error estimate have to be carried out to accomplish such a convergence
result. In our knowledge, this work will be the first to combine the following
three theoretical properties for a numerical scheme for the PNP system: (i)
unique solvability and positivity, (ii) energy stability, and (iii) optimal
rate convergence. A few numerical results are also presented in this article,
which demonstrates the robustness of the proposed numerical scheme.
