The self-consistent modeling of vacuum polarization due to virtual
electron-positron fluctuations is of relevance for many near term experiments
associated with high intensity radiation sources and represents a milestone in
describing scenarios of extreme energy density. We present a generalized
finite-difference time-domain solver that can incorporate the modifications to
Maxwell's equations due to vacuum polarization. Our multidimensional solver
reproduced in one dimensional configurations the results for which an analytic
treatment is possible, yielding vacuum harmonic generation and birefringence.
The solver has also been tested for two-dimensional scenarios where finite
laser beam spot sizes must be taken into account. We employ this solver to
explore different types of counter-propagating configurations that can be
relevant for future planned experiments aiming to detect quantum vacuum
dynamics at ultra-high electromagnetic field intensities
