  Discretizing a distribution function in a phase space for an efficient
quantum dynamics simulation is a non-trivial challenge, in particular for a
case that a system is further coupled to environmental degrees of freedom. Such
open quantum dynamics is described by a reduced equation of motion (REOM) most
notably by a quantum Fokker-Planck equation (QFPE) for a Wigner distribution
function (WDF). To develop a discretization scheme that is stable for numerical
simulations from the REOM approach, we employ a two-dimensional (2D)
periodically invariant system-bath (PISB) model with two heat baths. This model
is an ideal platform not only for a periodic system but also for a non-periodic
system confined by a potential. We then derive the numerically ''exact''
hierarchical equations of motion (HEOM) for a discrete WDF in terms of
periodically invariant operators in both coordinate and momentum spaces. The
obtained equations can treat non-Markovian heat-bath in a non-perturbative
manner at finite temperatures regardless of the mesh size. As demonstrations,
we numerically integrate the discrete QFPE for a 2D free rotor and harmonic
potential systems in a high-temperature Markovian case using a coarse mesh with
initial conditions that involve singularity.
