  In this work, we propose an efficient and accurate computational method to
evaluate the many-potential $\alpha\left(Z\alpha\right)^{n\ge3}$ vacuum
polarization density of hydrogen-like atoms within the finite-basis
approximation of the Dirac equation. To prove the performance of our
computational method, we choose to work with the one-electron
$_{\,\,\,92}^{238}\text{U}$ atom. In summary, we find that compliance with
charge conjugation symmetry is a priori required to obtain physical results
that are in line with our knowledge of the analytical problem. We also note
that the final numerical results are found to be in excellent agreement with
previous formal analytical (and numerical) evaluations that are limited to a
few simple nuclear distribution models. Our technique can be efficiently
implemented and evaluated in codes that solve the radial Dirac equation in the
finite basis set framework and allows the use of arbitrary (radial) nuclear
charge distribution. The obtained numerical results of the non-perturbative
vacuum polarization density automatically account for the extended nuclear size
effect. This method is hence of special importance for atomic Dirac problems
whose analytical Green's functions expressions are not at hand or have
relatively complicated analytical forms. Furthermore, we propose a vacuum
polarization density formula that forces compliance with charge conjugation
symmetry and can be used in cases where the relativistic basis violates this
symmetry, as is the case in most relativistic basis set programs. In addition,
we have shown that vector components of the vacuum polarization four-current
vanish in the case where the Dirac Hamiltonian is symmetric under time-reversal
symmetry.
