  Two finite difference schemes for the numerical treatment of the von Neumann
equation for the (2+1)D Dirac Hamiltonian are presented. Both utilize a
single-cone staggered space-time grid which ensures a single-cone energy
dispersion to formulate a numerical treatment of the mixed-state dynamics
within the von Neumann equation. The first scheme executes the time-derivative
according to the product rule for "bra" and "ket" indices of the density
operator. It therefore directly inherits all the favorable properties of the
difference scheme for the pure-state Dirac equation and conserves positivity.
The second scheme proposed here performs the time-derivative in one sweep. This
direct scheme is investigated regarding stability and convergence. Both schemes
are tested numerically for elementary simulations using parameters which
pertain to topological insulator surface states. Application of the schemes to
a Dirac Lindblad equation and real-space-time Green's function formulations are
discussed.
