  Nonequilibrium dynamics of highly-controlled quantum systems is a challenging
issue in statistical physics and quantum many-body physics, relevant to recent
experimental developments of analog and digital quantum simulations. In this
work, we develop a discrete phase-space approach for general SU($N$) spin
systems that utilizes cluster mean field equations, which capture non-trivial
quantum correlations inside each cluster, beyond the capability of the standard
discrete truncated Wigner approximation for individual classical spins. Our
formalism, based on a cluster phase-point operator, makes it possible to
realize scalable numerical samplings of cluster phase-space variables, where
the total number of noise variables for a direct product state is independent
of the choice of the separation into finite regions of clusters. We numerically
demonstrate that the cluster discrete truncated Wigner approximation (C-dTWA)
method can reproduce key results in a recent experiment on the correlation
propagation dynamics in a two dimensional Bose-Hubbard system. We also compare
the results of C-dTWA for clusters of $2\times 2$ sites with those of a
two-dimensional tensor network method and discuss that both approaches agree
very well in a short time region, where the energy is well conserved in the
tensor network simulations. Since we formulate the C-dTWA method in a general
form, it can be potentially applied to various dynamical problems in isolated
and open quantum systems even in higher dimensions.
