# Hybrid discrete-continuous truncated Wigner approximation
for driven, dissipative spin systems

Christopher D. Mink  Department of Physics and Research Center OPTIMAS,
University of Kaiserslautern, D-67663 Kaiserslautern, Germany  David Petrosyan
Institute of Electronic Structure and Laser, FORTH, GR-70013 Heraklion, Crete,
Greece  Michael Fleischhauer  Department of Physics and Research Center
OPTIMAS, University of Kaiserslautern, D-67663 Kaiserslautern, Germany

###### Abstract

We present a systematic approach for the semiclassical treatment of many-body
dynamics of interacting, open spin systems. Our approach overcomes some of the
shortcomings of the recently developed discrete truncated Wigner approximation
(DTWA) based on Monte Carlo sampling in a discrete phase space that improves
the classical treatment by accounting for lowest-order quantum fluctuations.
We provide a rigorous derivation of the DTWA by embedding it in a continuous
phase space, thereby introducing a hybrid discrete-continuous truncated Wigner
approximation. We derive a set of operator-differential mappings that yield an
exact equation of motion (EOM) for the continuous SU(2) Wigner function of
spins. The standard DTWA is then recovered by a systematic neglection of
specific terms in this exact EOM. The hybrid approach permits us to determine
the validity conditions and to gain a detailed understanding of the quality of
the approximation, paving the way for systematic improvements. Furthermore, we
show that the continuous embedding allows for a straightforward extension of
the method to open spin systems subject to dephasing, losses, and incoherent
drive, while preserving the key advantages of the discrete approach, such as a
positive definite Wigner distribution of typical initial states. We derive
exact stochastic differential equations for processes which cannot be
described by the standard DTWA due to the presence of non-classical noise. We
illustrate our approach by applying it to the dissipative dynamics of Rydberg
excitation of one-dimensional arrays of laser-driven atoms and compare it to
exact results for small systems.

##  I Introduction

The many-body dynamics of dissipative quantum spin systems is of key
importance in many areas of physics and technology. Its exact numerical
treatment is however extremely challenging, being restricted either to small
systems where the time evolution of the full many-body density matrix can be
simulated or to the classical limit of strong dephasing, which can be tackled
by Monte Carlo methods [1, 2]. Various approximation techniques have been
developed in the past, ranging from mean-field, cluster-mean-field [3], and
variational approaches [4], as well as field-theoretical descriptions within
the Keldysh formalism [5] to those based on matrix-product state (MPS)
expansions of the density matrix [6, 7] or variational MPS techniques [8].
More recently, a semi-classical approach based on Monte Carlo sampling of
spin-$\frac{1}{2}$ density matrices in a discrete phase space [9] – the
discrete truncated Wigner approximation (DTWA) [10] – has been developed. In
the DTWA, the many-body dynamics of spin-$\frac{1}{2}$ systems is described by
a set of classical equations of motion of the Cartesian spin components, while
interactions are accounted for at the mean-field level. Similar to the
truncated Wigner approach to Bose fields [11, 12, 13, 14], quantum
fluctuations are partially incorporated by sampling the spin components from
an initial quasi-probability distribution of a discrete phase space [9].

It was shown in Ref. [10] that the DTWA can reproduce rather accurately
collective observables and correlations of interacting spin-$\frac{1}{2}$
particles on short time scales and thus can be successfully applied to
problems of spin-squeezing [15, 16] and quantum quenches [17, 18, 19, 20]. The
agreement with exact solutions is particularly good for long-range
interactions, i.e., when each spin interacts with a large number of other
spins with comparable strength. While qualitatively such a behavior does not
come unexpected, since increasing the coordination number of interactions for
each spin improves the accuracy of the mean-field approximation, detailed
quantitative understanding of the applicability and limitations of the DTWA is
still missing. Yet, such an understanding is important in order to estimate
the quality of the DTWA when applied to larger systems, where exact benchmarks
can no longer be performed, or for the development of systematic improvements
[21].

Furthermore, while phase space methods for bosonic fields [11, 13, 12] can
easily incorporate dephasing or losses, so far there has been no first-
principles derivation of the DTWA for general open systems. Only very recently
a first attempt has been made to include dephasing and incoherent decay using
a phenomenological approach. In Ref. [22] dephasing was incorporated by
Markovian classical noise fields coupled to the $x,y$ components of the spins,
turning the deterministic equations of motion into stochastic ones with
multiplicative white noise. In contrast, decay and incoherent drive are
associated with non-classical noise in standard DTWA, which prevents numerical
simulations by stochastic differential equations. Further approximations were
therefore made [22, 23]. In an alternative approach of treating pure decay,
the authors of Ref. [24] proposed to first unravel the quantum master equation
and then apply a semiclassical approach. This open-system version of the DTWA
requires an ad hoc introduction of the spin norm as an additional degree of
freedom. For a single spin, very good agreement with exact results was found,
but the range of validity of the approach as discussed in Appendix D and its
applicability to systems with simultaneous dephasing and loss or to more
general reservoir couplings remains unclear.

Here we present an alternative, rigorous derivation of the DTWA based on an
embedding of the discrete Wigner representation in the continuous SU(2) phase
space. This hybrid discrete-continuous approach, which we term as DCTWA
(discrete-continuous truncated Wigner approximation), alleviates the
shortcomings of the standard DTWA, namely, it permits us to assess and
systematically improve the quality of the truncation approximation and allows
to extend the approach to treat both dephasing and decay. To this end, we
derive a set of operator-differential mappings which yield an exact equation
of motion (EOM) for the continuous SU(2) Wigner function. The truncation
approximation, which forms the basis of the DTWA, can then be associated with
a systematic neglection of higher order derivatives in this EOM. A scaling
analysis of the higher-order terms provides clear insight into the range of
validity of the approximation and paves the way for systematic improvements.
Furthermore, the continuous approach provides a straightforward extension to
open-system dynamics. In sharp contrast to the DTWA, the diffusion matrix
associated with spin decay is always positive definite for the continuous
Wigner distribution. Yet, typical initial states have continuous Wigner
distributions which are not positive definite and thus cannot be treated as
probability distributions, preventing a Monte Carlo sampling. We here show
that exploiting the connection between discrete and continuous
representations, together with the gauge degree of freedom of the continuous
Wigner function of spins, allows for both a Monte Carlo sampling of initial
states as well as a simulation of time evolution by stochastic differential
equations. Moreover, just like the approach of Ref. [24] the dynamics of a
single driven-dissipative spin is treated exactly (see Appendix D).

The paper is organized as follows. In Sec. II we briefly review Wootters’
discrete Wigner function formalism for spins [9] and the standard DTWA based
on that formalism [10]. In Sec. III we summarize the continuous SU(2) Wigner
representation of spins [25, 26]. In Sec. IV we derive a mapping between the
discrete and continuous representations and develop a hybrid discrete-
continuous version, where the sampling of the initial distribution encoding
quantum fluctuations to lowest order is performed on a discrete space, while
the time evolution is performed in the continuous representation. In Sec. V we
illustrate the performance of the DCTWA by applying it to the dissipative
dynamics of Rydberg excitation of one-dimensional arrays of laser-driven atoms
[27, 28] and compare it to exact results for small atomic systems. Our
findings are summarized in Sec. VI.

##  II Discrete Wigner function for spins and Truncation approximation

###  II.1 Quantum systems with continuous degrees of freedom

The density operator $\hat{\rho}$ of a quantum system with continuous degrees
of freedom can equivalently be represented by the Wigner function. For
example, a quantum particle with position $\hat{q}$ and momentum $\hat{p}$ has
a corresponding quasiprobability distribution $W(q,p)$ of continuous scalars
$q,p$. The connection between the Hilbert space and Wigner phase space is
given by

$\hat{\rho}=\iint dq\,dp\,\,W(q,p)\,\hat{A}(q,p),$ (1)

where

$\hat{A}(q,p)=\int\\!dy\,\Bigl{|}q-\frac{y}{2}\Bigr{\rangle}\Bigl{\langle}q+\frac{y}{2}\Bigr{|}\,e^{-i\frac{py}{\hbar}}$

are continuous phase-point operators [29, 30, 31]. Expectation values of any
observable $\hat{O}$ can then be calculated by taking the statistical average
of the corresponding Weyl symbol $\mathcal{O}$ with respect to the Wigner
function

$\langle\hat{O}(\hat{q},\hat{p})\rangle=\iint dq\,dp\,\,W(q,p)\ {\cal
O}(p,q),$ (2)

over the whole phase space, where the Weyl symbol ${\cal O}$ is defined via

$\displaystyle{\cal O}(p,q)$ $\displaystyle=$
$\displaystyle\text{Tr}\left\\{\hat{A}(q,p)\hat{O}(\hat{q},\hat{p})\right\\}$
(3) $\displaystyle=$
$\displaystyle\int\\!dy\,\Bigl{\langle}q+\frac{y}{2}\Bigr{|}\hat{O}(\hat{q},\hat{p})\Bigl{|}q-\frac{y}{2}\Bigr{\rangle}\,e^{-i\frac{py}{\hbar}}.$

###  II.2 Wigner function for a two-state quantum system

The concept of phase space representations can be extended to quantum systems
with a finite-dimensional Hilbert space [9] such as spin-$\frac{1}{2}$
systems. To this end, Wootters [9] introduced four discrete phase-point
operators

$\displaystyle\hat{A}_{\bm{\alpha}}=$
$\displaystyle\frac{1}{2}\left(\hat{\mathbbm{1}}_{2}+\bm{r}_{\bm{\alpha}}\hat{\bm{\sigma}}\right),$
(4)

where
$\hat{\bm{\sigma}}=(\hat{\sigma}^{x},\hat{\sigma}^{y},\hat{\sigma}^{z})^{T}$
is the vector of Pauli matrices, $\bm{\alpha}=(\alpha_{1},\alpha_{2})$ with
elements $\alpha_{i}\in\\{0,1\\}$ and the vector $\bm{r}_{\bm{\alpha}}$ is
given by the discrete set of points

$\displaystyle\bm{r}_{\bm{\alpha}}=$
$\displaystyle\Bigl{(}(-1)^{\alpha_{2}},(-1)^{\alpha_{1}+\alpha_{2}},(-1)^{\alpha_{1}}\Bigr{)}$
(5)

on the sphere with radius $r=\sqrt{3}$. The phase-point operators have unit
trace $\text{Tr}(\hat{A}_{\bm{\alpha}})=1$, are orthogonal
$\frac{1}{2}\text{Tr}(\hat{A}_{\bm{\alpha}}\hat{A}_{\bm{\beta}})=\delta_{\bm{\alpha}\bm{\beta}}$,
and form a basis. The density operator $\hat{\rho}$ can then be expanded as

$\displaystyle\hat{\rho}=\sum_{\bm{\alpha}}W_{\bm{\alpha}}\hat{A}_{\bm{\alpha}},$
(6)

where the weights
$W_{\bm{\alpha}}=\frac{1}{2}\text{Tr}[\hat{A}_{\bm{\alpha}}\hat{\rho}]$ are
called the discrete Wigner function and are real by construction.

Note that since $r=\sqrt{3}>1$, the phase-point operators
$\hat{A}_{\bm{\alpha}}$ are not positive definite and therefore are not
density matrices themselves. As a consequence, the discrete Wigner function
$W_{\bm{\alpha}}$ is in general not positive, but is normalized
$\text{Tr}\,\hat{\rho}=\sum_{\bm{\alpha}}W_{\bm{\alpha}}=1$. Just as in the
continuous case, the Wigner function represents only a quasi-probability
distribution, therefore preventing a stochastic (Monte Carlo) sampling of
arbitrary states. But for certain relevant quantum states, all Wigner
coefficients are positive and can therefore be interpreted as proper
probabilities. For example, the fully-polarized spin states

$\displaystyle|\\!\uparrow\rangle\langle\uparrow\\!|$
$\displaystyle=\frac{1}{2}(\hat{A}_{00}+\hat{A}_{01})$
$\displaystyle\Longleftrightarrow W_{00}=W_{01}=\frac{1}{2},W_{10}=W_{11}=0,$
(7a) $\displaystyle|\\!\downarrow\rangle\langle\downarrow\\!|$
$\displaystyle=\frac{1}{2}(\hat{A}_{10}+\hat{A}_{11})$
$\displaystyle\Longleftrightarrow W_{00}=W_{01}=0,W_{10}=W_{11}=\frac{1}{2},$
(7b)

can be represented as equally weighted classical mixtures of $\hat{A}_{00}$
and $\hat{A}_{01}$ (for $|\\!\uparrow\rangle\langle\uparrow\\!|$), and
$\hat{A}_{10}$ and $\hat{A}_{11}$ (for
$|\\!\downarrow\rangle\langle\downarrow\\!|$).

###  II.3 Equations of motion

It follows from the von Neumann equation
$\frac{d}{dt}\hat{\rho}(t)=-i[\hat{H}_{0},\hat{\rho}(t)]$ ($\hbar=1$) for the
unitary time evolution under some single-particle Hamiltonian $\hat{H}_{0}$
that the evolution equations for the phase-point operators are

$\displaystyle\frac{d}{dt}\hat{A}_{\bm{\alpha}}(t)=-i[\hat{H}_{0},\hat{A}_{\bm{\alpha}}(t)].$
(8)

Note that these equations describe the time evolution in the Schrödinger
picture, as the phase point operators correspond to the quantum state. Since
each of the four phase-point operators $\hat{A}_{\bm{\alpha}}$ is a linear
combination of the Pauli spin matrices
$\hat{\sigma}^{x},\hat{\sigma}^{y},\hat{\sigma}^{z}$, their time evolution can
be expressed in terms of a time-dependent vector
$\bm{s}_{\bm{\alpha}}(t)=\Bigl{(}s^{x}_{\bm{\alpha}}(t),s^{y}_{\bm{\alpha}}(t),s^{z}_{\bm{\alpha}}(t)\Bigr{)}$,

$\displaystyle\hat{A}_{\bm{\alpha}}(t)=\frac{1}{2}\left(\hat{\mathbbm{1}}_{2}+\bm{s}_{\bm{\alpha}}(t)\hat{\bm{\sigma}}\right).$
(9)

Here $s^{\mu}(t)$ ($\mu=x,y,z$) obey the equations of motion for a classical
spin (see Appendix A),

$\displaystyle\dot{s}^{\mu}(t)=\bigl{\\{}s^{\mu},{\cal
H}_{0}\bigr{\\}}_{P}=2\sum_{\nu,\lambda}\epsilon_{\mu\nu\lambda}\,s^{\lambda}\,\frac{\partial{\cal
H}_{0}}{\partial s^{\nu}},$ (10)

where the classical Hamiltonian ${\cal H}_{0}({\bm{s}})$ is the Weyl symbol
corresponding to $\hat{H}_{0}$ and we have dropped the index $\bm{\alpha}$ for
notational simplicity. The initial conditions of the vector
$\bm{s}_{\bm{\alpha}}$ are: $\bm{s}_{\bm{\alpha}}(t=0)=\bm{r}_{\bm{\alpha}}$.
Evaluating the expectation value of any observable
$\langle\hat{O}(\hat{\bm{\sigma}}(t))\rangle$ then amounts to calculating the
time evolution of the Weyl symbol ${\cal
O}({\bm{s}}(t))=\textrm{Tr}\left\\{\hat{O}\hat{A}_{\bm{\alpha}}(t)\right\\}$
averaged over the initial discrete Wigner distribution

$\langle\hat{O}(\hat{\bm{\sigma}}(t))\rangle=\sum_{\bm{\alpha}}W_{\bm{\alpha}}(0)\,{\cal
O}({\bm{s}}(t)).$ (11)

If the initial state $\hat{\rho}(0)$ has positive Wigner coefficients
$W_{\bm{\alpha}}$, as for the state $|\negmedspace\uparrow\rangle$ or
$|\negmedspace\downarrow\rangle$ in Eqs. (7), Eq. (11) can be evaluated as a
classical average of the Weyl symbols weighted by the initial probabilities
$W_{\bm{\alpha}}(0)$. The ”quantumness” of the spin is then captured by the
averaging over the initial conditions. This is possible without further
approximation since we consider a single spin-$\frac{1}{2}$ particle without
interactions and decay, in which case the Heisenberg equations of motion are
also linear.

###  II.4 Interacting spin systems and DTWA

The utility of the discrete phase-space approach stems from its application to
many interacting spins, which in general have complicated and intractable
many-body dynamics. The phase space representation is the basis of the DTWA
[10], allowing it to take lowest-order quantum fluctuations into account and
thus going beyond the mean-field level. Below, we outline the essence of the
DTWA.

The density matrix for $N$ interacting spins can be faithfully represented by
$4^{N}$ discrete Wigner coefficients and the corresponding phase point
operators

$\displaystyle\hat{\rho}=\sum_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}W_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}\,\hat{A}_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}.$
(12)

As for a single spin, the dynamics of $\hat{\rho}$ can then be mapped onto the
dynamics of the $4^{N}$ phase point operators
$\hat{A}_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}(t)$, whose general evaluation
in time requires exponentially increasing effort in $N$. The key approximation
in the DTWA that makes the many-body problem tractable is the factorization of
the phase point operators at all times:

$\hat{A}_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}(t)\approx\hat{A}_{1,\bm{\alpha_{1}}}(t)\otimes\hat{A}_{2,\bm{\alpha_{2}}}(t)\otimes\ldots\otimes\hat{A}_{N,\bm{\alpha_{N}}}(t).$
(13)

The many-body density operator at time $t$ is then explicitly given by

$\displaystyle\hat{\rho}(t)=\sum_{\bm{\alpha}_{1},\dots\bm{\alpha}_{N}}\hat{A}_{1,\bm{\alpha}_{1}}(t)\otimes\ldots\otimes\hat{A}_{N,\bm{\alpha}_{N}}(t)W_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}(0),$
(14)

which, for a factorized initial state, is further simplified to
$W_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}(0)=\prod_{n}W_{\bm{\alpha}_{n}}(0)$.

Substituting the factorization ansatz of Eq. (13) into the equation of motion
with a many-body Hamiltonian $\hat{H}$,

$\displaystyle\frac{d}{dt}\hat{A}_{\bm{\alpha}_{1},\dots\bm{\alpha}_{N}}(t)=-i[\hat{H},\hat{A}_{\bm{\alpha}_{1},\dots\bm{\alpha}_{N}}(t)],$
(15)

and using the orthogonality of the phase point operators, one finds that the
time evolution of the single-spin operators is governed by

$\displaystyle\frac{d}{dt}\hat{A}_{j,\bm{\alpha}_{j}}(t)$ $\displaystyle=$
$\displaystyle-i\bigl{[}\hat{H}_{j}^{\textrm{MF}},\hat{A}_{j,\bm{\alpha}_{j}}\bigr{]}.$
(16)

Here the mean-field Hamiltonian $\hat{H}_{j}^{\textrm{MF}}$ for $j$th spin is
obtained by replacing all the other spin operators $\hat{\bm{\sigma}}_{l}$ by
the corresponding c-numbers $\bm{s}_{l}^{\bm{\alpha}}(t)$. Using the
decomposition of the phase point operators in Eq. (9), we thus find that
$s_{j}^{\mu}(t)$ $(\mu=x,y,z)$ obey the classical equations of motion with the
classical many-body Hamiltonian,

$\displaystyle\dot{s}_{j}^{\mu}(t)=\bigl{\\{}s_{j}^{\mu},{\cal
H}\bigr{\\}}_{P}=2\sum_{\nu}^{\lambda}\epsilon_{\mu\nu\lambda}\,s_{j}^{\lambda}\,\frac{\partial{\cal
H}}{\partial s_{j}^{\nu}},$ (17)

with initial conditions
$\bm{s}_{j,\bm{\alpha}_{j}}(t=0)=\bm{r_{j,\bm{\alpha}_{j}}}$. As in the case
of a single spin, all quantum mechanical observables are obtained from the
solutions of the classical (mean-field) equations (17) and averaging over the
initial Wigner distribution

$\displaystyle\langle\hat{O}(\hat{\bm{\sigma}}_{1},$
$\displaystyle\dots,\hat{\bm{\sigma}}_{N};t)\rangle=$
$\displaystyle\sum_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}\\!\\!\\!W_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}(0)\,{\cal
O}(\bm{s}_{1}(t),\dots,\bm{s}_{N}(t)).$ (18)

It is worth noting that, despite the mean-field character of the DTWA, the
quantum nature of the problem is taken into account to some extent in two
ways. Firstly, as in the case of a single spin, the averaging over the initial
conditions for every spin captures lowest-order quantum fluctuations on the
single-particle level. Secondly, information about initial (quantum)
correlations between the spins is, in principle, contained in
$W_{\bm{\alpha}_{1},\dots,\bm{\alpha}_{N}}(0)$. Hence, the DTWA goes beyond a
mean-field description of the many-spin problem, as we will illustrate below
with a specific example of an optically driven lattice of Rydberg atoms. This
and the fact that it is computationally not much more expensive than a mean-
field calculation makes the DTWA an appealing semiclassical approach. However,
the derivation of the DTWA is based on the heuristic approximation of
factorizing the many-body phase point operators. The quantitative
characterization and range of validity of this approximation are not well
defined, and the level of improvement over a simple mean-field approach is not
clear.

###  II.5 Open spin systems and DTWA

It is tempting to extend the above procedure to open spin systems coupled to
Markovian reservoirs. This is described by the Lindblad master equation

$\displaystyle\frac{d}{dt}\hat{\rho}=-i\left[\hat{H},\hat{\rho}\right]+\frac{1}{2}\sum_{\mu}\left(2\hat{L}_{\mu}\hat{\rho}\hat{L}_{\mu}^{\dagger}-\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}\hat{\rho}-\hat{\rho}\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}\right),$
(19)

where the $\hat{L}_{\mu}$ are the Lindblad generators. Substituting the
representation of $\hat{\rho}(t)$ in terms of the discrete phase-point
operators $\hat{A}_{\bm{\alpha}}(t)$, one recognizes however that, in general,
they do not remain orthogonal under non-unitary time evolution. This prevents
a straightforward extension of the DTWA to open systems and one has to resort
to approaches where the Lindblad master equation is effectively generated,
e.g., by coupling to classical noise fields [22] or by an unraveling procedure
[24].

##  III Continuous Wigner function for spins

One of the key advantages of continuous phase-space representations of quantum
states is the connection to stochastic processes via equations of motion which
are partial differential equations and — with additional approximations — are
of Fokker-Planck type [32]. This permits efficient numerical simulation of the
system dynamics in terms of stochastic differential equations. Despite its
successful application in quantum optics, continuous Wigner functions had not
been widely applied to finite dimensional quantum systems, such as spins.
Following the approach of Ref. [25], Tilma et. al. [26] defined a Wigner
distribution $W_{\hat{\rho}}(\Omega)$ for a quantum state $\hat{\rho}$ over a
continuous phase space characterized by parameters $\Omega$, provided there
exists a kernel $\hat{A}(\Omega)$ that generates $W_{\hat{\rho}}(\Omega)$
according to the generalized Weyl rule
$W_{\hat{\rho}}(\Omega)=\text{Tr}[\hat{\rho}\hat{A}(\Omega)]$ and that
satisfies the Stratonovich-Weyl correspondence [33].

The kernel operator $\hat{A}(\Omega)$ and the set of coordinates $\Omega$ are
not unique. For a spin-$\frac{1}{2}$ system, there exists, in particular, a
representation of the state $\hat{\rho}$ through

$\displaystyle\hat{A}(\theta,\phi)$ $\displaystyle=$ $\displaystyle
U(\theta,\phi,\psi)\hat{A}_{0}U^{\dagger}(\theta,\phi,\psi),$ (20)

where
$\hat{A}_{0}=\frac{1}{2}(\hat{\mathbbm{1}}_{2}-\sqrt{3}\hat{\sigma}^{z})$, and
$U(\theta,\phi,\psi)=e^{i\hat{\sigma}^{z}\phi/2}e^{i\hat{\sigma}^{y}\theta/2}e^{i\hat{\sigma}^{z}\psi/2}$,
are the SU(2) rotation operators with the Euler angles $(\theta,\phi,\psi)$
that span the continuous phase space of the surface of a sphere,
$\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$. Note that in Eq. (20) the dependence
on the angle $\psi$ drops out and one finds

$\displaystyle\hat{A}(\theta,\phi)=$
$\displaystyle\frac{1}{2}[\hat{\mathbbm{1}}_{2}+\bm{s}(\theta,\phi)\hat{\bm{\sigma}}]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}1-\sqrt{3}\cos\theta&&\sqrt{3}e^{i\phi}\sin\theta\\\
\sqrt{3}e^{-i\phi}\sin\theta&&1+\sqrt{3}\cos\theta\end{pmatrix},$ (21)
$\displaystyle=$
$\displaystyle\frac{\sqrt{4\pi}}{2}\Bigl{(}\text{Y}_{00}(\theta,\phi)\,\hat{\mathbbm{1}}_{2}-\text{Y}_{10}(\theta,\phi)\,\hat{\sigma}^{z}\Bigr{)}$
$\displaystyle+\sqrt{2\pi}\frac{\text{Y}_{1,-1}(\theta,\phi)-\text{Y}_{11}(\theta,\phi)}{2}\hat{\sigma}^{x}$
$\displaystyle-i\sqrt{2\pi}\frac{\text{Y}_{1,-1}(\theta,\phi)+\text{Y}_{11}(\theta,\phi)}{2}\hat{\sigma}^{y}.$

Here the c-number vector

$\displaystyle\bm{s}(\theta,\phi)=\sqrt{3}(\sin\theta\cos\phi,-\sin\theta\sin\phi,-\cos\theta)^{T}$
(22)

is a representation of the surface of the sphere with radius $\sqrt{3}$ and
the $\text{Y}_{lm}(\theta,\phi)$ are the spherical harmonics.

The relation between any observable $\hat{O}(\hat{\bm{\sigma}})$ for a
spin-$\frac{1}{2}$ system and the corresponding Weyl symbol follows from the
simple algebra of the Pauli matrices. Since any function $\hat{O}$ of spin
operators can be written as a linear superposition of the unity matrix and the
Pauli matrices $\hat{\sigma}^{\mu}$ ($\mu=(x,y,z)$),

$\hat{O}(\hat{\bm{\sigma}})=a_{0}\hat{\mathbbm{1}}_{2}+a_{y}\hat{\sigma}^{x}+a_{y}\hat{\sigma}^{y}+a_{z}\hat{\sigma}^{z},$
(23)

its Weyl symbol is obtained by replacing the spin operators by the
corresponding components of the classical spin vector $\bm{s}(\theta,\phi)$,

${\cal
O}\bigl{(}\bm{s}(\theta,\phi)\bigr{)}=\text{Tr}\left\\{\hat{A}(\theta,\phi)\hat{O}(\hat{\bm{\sigma}})\right\\}=O({\bm{s}}).$
(24)

An arbitrary spin state $\hat{\rho}$ can now be expressed via the continuous
Wigner function $W(\theta,\phi)\in\mathbb{R}$ and vice versa as

$\displaystyle\hat{\rho}=\int d\Omega\,\hat{A}(\theta,\phi)W(\theta,\phi),$
(25a) $\displaystyle
W(\theta,\phi)=\text{Tr}[\hat{A}(\theta,\phi)\hat{\rho}],$ (25b)

where $\int d\Omega=\int_{0}^{\pi}d\theta\sin\theta\int_{0}^{2\pi}d\phi/2\pi$
and $\hat{A}=\hat{A}^{\dagger}$. The Wigner function is normalized $\int
d\Omega\,W(\theta,\phi)=\text{Tr}\,\hat{\rho}=1$, but is in general not
positive, i.e., is a quasi-probability distribution. Note that $\int
d\Omega\,\hat{A}(\theta,\phi)=\hat{\mathbbm{1}}_{2}$.

The continuous Wigner representations for fully-polarized spin states are

$\displaystyle W_{\uparrow}(\theta,\phi)=$
$\displaystyle\text{Tr}[\hat{A}(\theta,\phi)|\\!\uparrow\rangle\langle\uparrow\\!|]=\frac{1}{2}(1-\sqrt{3}\cos\theta),$
(26a) $\displaystyle W_{\downarrow}(\theta,\phi)=$
$\displaystyle\text{Tr}[\hat{A}(\theta,\phi)|\\!\downarrow\rangle\langle\downarrow\\!|]=\frac{1}{2}(1+\sqrt{3}\cos\theta).$
(26b)

These Wigner functions are not positive semi-definite and therefore cannot be
approximated faithfully by a Monte Carlo sampling. In fact all pure spin
states have non-positive continuous Wigner functions and cannot be sampled.
Note, moreover, that the surface element $d\Omega=\sin\theta d\theta
d\phi/2\pi$ is nonlinear in $\theta$, which prevents straightforward
derivation of a Fokker-Planck type equation of motion for $W(\theta,\phi)$.
Both problems can be resolved, however, by a hybrid discrete-continuous
approach as outlined in the following section.

##  IV Hybrid discrete-continuous truncated Wigner approximation

The shortcomings of the standard DTWA can be overcome by the hybrid DCTWA
approach in which the sampling over initial conditions is performed for the
Cartesian spin components as in the standard DTWA, while the time evolution is
performed in an angular representation using a ”flattened” continuous Wigner
distribution that leads to a Fokker-Planck type EOM.

###  IV.1 Time evolution of continuous Wigner function for a single spin

We first derive a mapping from the Lindblad master equation (19) of the
density operator $\hat{\rho}$ for a single spin to an equation of motion of
the continuous Wigner function $W(\theta,\phi)$. To this end, we employ the
decomposition (25) of $\hat{\rho}$ into continuous phase point operators and
note that
$\hat{A}(\theta,\phi),\frac{\partial}{\partial\theta}\hat{A}(\theta,\phi),\frac{\partial}{\partial\phi}\hat{A}(\theta,\phi)$
and $\frac{\partial^{2}}{\partial\phi^{2}}\hat{A}(\theta,\phi)$ are linearly
independent and thus form a basis for all $2\times 2$ matrices. We can,
therefore, express any operator acting on $\hat{\rho}$ by its action on a
phase-point operator using an operator-differential identity, such as

$\displaystyle\hat{\sigma}^{z}\hat{A}(\theta,\phi)=$
$\displaystyle\bigg{[}-\sqrt{3}\cos\theta+\frac{3\sin\theta-2\csc\theta}{\sqrt{3}}\frac{\partial}{\partial\theta}$
$\displaystyle-i\frac{\partial}{\partial\phi}-\frac{2\cot\theta\csc\theta}{\sqrt{3}}\frac{\partial^{2}}{\partial\phi^{2}}\bigg{]}\hat{A}(\theta,\phi),$
(27)

which can be proven by a straightforward evaluation of both sides.

Following the standard procedure of phase-space approaches, we insert these
mappings into the master equation (19), which, after partial integration,
leads to a dynamical equation for $W(\theta,\phi)$. In contrast to typical
quantum optical problems of interacting bosonic fields, however, the EOM is
not of (generalized) Fokker-Planck type, due to the nonlinear integration
measure on the sphere $d\Omega=\sin\theta\,d\theta d\phi/2\pi$. To circumvent
this problem, we introduce the ”flattened” Wigner function (FWF)

$\chi(\theta,\phi)=\frac{\sin\theta}{2\pi}\,W(\theta,\phi),$ (28)

defined on a stripe $\theta\in[0,\pi]$, $\phi\in[0,2\pi)$, which is periodic
in $\phi$ with period $2\pi$. The distinction between the Wigner function and
the FWF is critical for the time evolution, since the EOM of $W(\theta,\phi)$
contains in general terms that cannot be expressed through partial
derivatives. In contrast, the time evolution of the FWF can be expressed as a
Fokker-Planck equation (FPE) (provided higher than second-order derivatives
can be neglected), which in turn allows for an efficient numerical evaluation
of the system dynamics by stochastic differential equations.

More formally, this reinterpretation can be avoided by introducing the
contravariant coordinate vector $(x^{1},x^{2})=(\theta,\phi)$ and metric
tensor $g_{\mu\nu}$ of the curved, i.e., spherical, phase space which is given
by

$\displaystyle g=\frac{1}{2\pi}\begin{pmatrix}1&0\\\
0&\sin^{2}\theta\end{pmatrix}.$ (29)

The infinitesimal volume element is generated by
$\sqrt{\det(g)}=\sin\theta/2\pi$. If we were to use covariant derivatives

$\displaystyle\nabla_{\mu}=\frac{1}{\sqrt{\det(g)}}\frac{\partial}{\partial
x^{\mu}}\sqrt{\det(g)}$ (30)

instead of plain derivatives, the introduction of the FWF would not be
necessary and instead diffusion processes on the spherical surface would be
obtained. Since both approaches yield the same differential equations and
therefore identical results, we choose the one based on plain derivatives and
the FWF.

As opposed to the derivation in Sec. II, we now include the time dependence in
the Wigner function, or the flattened Wigner function, while keeping the
phase-point operators constant,

$\displaystyle\hat{\rho}(t)$ $\displaystyle=$ $\displaystyle\iint
d\theta\,d\phi\,\chi(\theta,\phi;t)\,\hat{A}(\theta,\phi).$ (31)

We substitute this expression into the Lindblad equation, evaluate the action
of spin operators on the kernel operators as in Eqs. (27) and (66), and
integrate by parts in order to let the derivatives act on the FWF instead of
on the kernel. Since the $(\theta,\phi)$ space is compact, the surface terms
vanish and we obtain

$\displaystyle\hat{\sigma}^{z}\hat{\rho}=$ $\displaystyle\int d\theta
d\phi\,\hat{A}(\theta,\phi)\bigg{[}-\sqrt{3}\cos\theta-\frac{\partial}{\partial\theta}\frac{3\sin\theta-2\csc\theta}{\sqrt{3}}$
$\displaystyle+\frac{\partial}{\partial\phi}i-\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\cot\theta\csc\theta}{\sqrt{3}}\bigg{]}\chi(\theta,\phi).$
(32)

This can be understood as a mapping between Hilbert space and phase space. We
abbreviate this relation to

$\displaystyle\hat{\sigma}^{z}\hat{\rho}\leftrightarrow$
$\displaystyle\bigg{[}-\sqrt{3}\cos\theta-\frac{\partial}{\partial\theta}\frac{3\sin\theta-2\csc\theta}{\sqrt{3}}$
$\displaystyle+\frac{\partial}{\partial\phi}i-\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\cot\theta\csc\theta}{\sqrt{3}}\bigg{]}\chi(\theta,\phi).$
(33)

Note that if we let the derivatives act on the Wigner function instead of on
the FWF, the contributions from the factor $\sin\theta$ would produce
different and far more involved mappings. A complete list of mappings for all
spin operators is given in Appendix B.

The mappings allow us to investigate the dynamics of a given system in terms
of the Wigner phase space. As an example, the unitary dynamics of a single
spin is fully governed by the differential contributions

$\displaystyle-i[\hat{\sigma}^{x},\hat{\rho}]\,\leftrightarrow$
$\displaystyle\,2\left(\frac{\partial}{\partial\theta}\sin\phi+\frac{\partial}{\partial\phi}\cot\theta\cos\phi\right)\chi,$
(34a) $\displaystyle-i[\hat{\sigma}^{y},\hat{\rho}]\,\leftrightarrow$
$\displaystyle\,2\left(\frac{\partial}{\partial\theta}\cos\phi-\frac{\partial}{\partial\phi}\cot\theta\sin\phi\right)\chi,$
(34b) $\displaystyle-i[\hat{\sigma}^{z},\hat{\rho}]\,\leftrightarrow$
$\displaystyle\,2\frac{\partial}{\partial\phi}\chi.$ (34c)

If we additionally consider the Lindblad master equation with a set of
Lindblad generators $\hat{L}_{\mu}$ that describe couplings to Markovian
reservoirs, we obtain an equation of the general form

$\displaystyle\int d\theta d\phi\,\hat{A}(\theta,\phi)\frac{\partial}{\partial
t}\chi=$ $\displaystyle\int d\theta
d\phi\,\hat{A}(\theta,\phi)\mathcal{L}\chi,$ (35)

where $\mathcal{L}$ is a differential operator. We thus see that a flattened
Wigner function $\chi(\theta,\phi,t)$ satisfying the partial differential
equation (PDE)

$\displaystyle\frac{\partial}{\partial t}\chi(\theta,\phi,t)=$
$\displaystyle\mathcal{L}\chi(\theta,\phi,t)$ (36)

faithfully represents the density matrix of Eq. (31).

Up to this point, no approximation has been made and the PDE describes the
full quantum problem. One type of PDE that commonly occurs in the context of
single-spin dynamics is the Fokker-Planck equation given by

$\displaystyle\frac{\partial}{\partial t}\chi(\bm{x},t)=$
$\displaystyle-\sum_{n}\frac{\partial}{\partial
x_{n}}\left[A_{n}(\bm{x},t)\chi(\bm{x},t)\right]$
$\displaystyle+\sum_{mn}\frac{\partial^{2}}{\partial x_{m}\partial
x_{n}}\left[D_{mn}(\bm{x},t)\chi(\bm{x},t)\right],$ (37)

where $\bm{x}=(\theta,\phi)^{T}$. We call $\bm{A}(\bm{x},t)$ the drift vector
and $\bm{D}(\bm{x},t)=\frac{1}{2}\bm{B}(\bm{x},t)^{T}\bm{B}(\bm{x},t)$ the
(positive-semidefinite) diffusion matrix. Instead of solving the Fokker-Planck
equation, one can then solve the corresponding Itô stochastic differential
equation (SDE)

$\displaystyle d\bm{x}(t)=$ $\displaystyle\bm{A}(t)dt+\bm{B}(t)d\bm{W},$ (38)

where $d\bm{W}=(dW_{\theta},dW_{\phi})^{T}$ is a multivariate differential
Wiener process [34].

This is useful for a numerical integration as one can sample a sufficient
number of initial states $\bm{x}(t=0)$ from a given positive semidefinite
Wigner function, determine their time evolution according to Eq. (38) and
calculate averages with respect to all trajectories to obtain a desired
observable:

$\displaystyle\left\langle\hat{O}(t)\right\rangle=$ $\displaystyle\int d\theta
d\phi\,\chi(\theta,\phi,t)\text{Tr}\left[\hat{O}\hat{A}(\theta,\phi)\right]$
$\displaystyle=$
$\displaystyle\overline{\biggl{.}\text{Tr}\left[\hat{O}\hat{A}(\theta(t),\phi(t))\right]}$
(39) $\displaystyle\approx$
$\displaystyle\frac{1}{M}\sum_{m=1}^{M}\text{Tr}\left[\hat{O}\hat{A}(\theta_{m}(t),\phi_{m}(t))\right],$

where $\theta_{m}(t),\phi_{m}(t)$ are the time-evolved trajectories and
$M<\infty$ is the number of evolved trajectories and the overline denotes
averaging over trajectories. As an example, in Fig. 1 we show two stochastic
trajectories for the evaluation of the the dynamics of a single
spin-$\frac{1}{2}$.

Figure 1: The phase space of a single spin-$\frac{1}{2}$, or a qubit, is
defined via the parametrization of the phase-point operator
$\hat{A}(\theta,\phi)$ in terms of angles $(\theta,\phi)$ on a sphere with
radius $\sqrt{3}$. The time evolution of a state $\hat{\rho}$ in Hilbert space
is equivalently expressed by an ensemble average of stochastic trajectories
$\\{\theta_{i}(t),\phi_{i}(t)\\}$. Here, two exemplary stochastic paths (blue
and orange lines), starting from the same discrete initial condition, are
shown for the decay process of Eqs. (48).

###  IV.2 Many spins and truncated Wigner approximation

The above results for a single spin can be generalized to a system of $N>1$
spins by constructing the phase-space representation of the many-body density
operator

$\displaystyle\hat{\rho}(t)=$ $\displaystyle\left[\prod_{n=1}^{N}\int
d\Omega_{n}\hat{A}(\theta_{n},\phi_{n})\right]W(t,\bm{\theta},\bm{\phi}),$
(40)

where $\bm{\theta}=(\theta_{1},\ldots,\theta_{N})^{T}$,
$\bm{\phi}=(\phi_{1},\ldots,\phi_{N})^{T}$ and
$d\Omega_{n}=d\theta_{n}\sin\theta_{n}d\phi_{n}/2\pi$ pertains to the $n$th
spin. A Pauli matrix $\hat{\sigma}_{n}^{\xi}$ ($\xi=x,y,z$) in the subspace of
the $n$th spin acting on
$\hat{A}(\bm{\theta},\bm{\phi})=\prod_{n}\hat{A}(\theta_{n},\phi_{n})$
generates the same operator-differential identities as in Eqs. (66). The FWF
then becomes
$\chi(\bm{\theta},\bm{\phi})=\left(\prod_{n=1}^{N}\frac{\sin\theta_{n}}{2\pi}\right)W(\bm{\theta},\bm{\phi})$.

Let us now consider interactions between the spins, e.g. of Ising type:

$\displaystyle\hat{H}=-\frac{1}{2}\sum_{m<n}J_{mn}\hat{\sigma}_{m}^{z}\hat{\sigma}_{n}^{z}.$
(41)

In analogy to the single spin approach, we derive equations for the time
evolution of the phase-point operators first. The interactions then generate
mixed-spin derivatives,

$\displaystyle\frac{d}{dt}\hat{A}(\bm{\theta},\bm{\phi})=$
$\displaystyle-\frac{i}{2}\sum_{m<n}J_{mn}[\hat{\sigma}_{m}^{z}\hat{\sigma}_{n}^{z},\hat{A}(\bm{\theta},\bm{\phi})]$
$\displaystyle=$ $\displaystyle\sum_{n}\sum_{m\neq
n}J_{mn}\bigg{(}\sqrt{3}\cos\theta_{m}\frac{\partial}{\partial\phi_{n}}+\frac{2\csc\theta_{m}-3\sin\theta_{m}}{\sqrt{3}}\frac{\partial^{2}}{\partial\theta_{m}\partial\phi_{n}}+\frac{2\cot\theta_{m}\csc\theta_{m}}{\sqrt{3}}\frac{\partial^{3}}{\partial\phi_{m}^{2}\partial\phi_{n}}\bigg{)}\hat{A}(\bm{\theta},\bm{\phi}).$
(42)

If only terms with derivatives in the same spin were present, the time
evolution of an initially factorized operator
$\hat{A}(t=0)=\prod_{n}\hat{A}_{n}(t=0)$ would preserve the factorization. But
the terms with mixed-spin derivatives, such as the second and third term of
Eq. (42) entangle the spins. Hence, in order to reproduce the DTWA, we have to
discard all terms with mixed spin derivatives, $n\neq m$. We note that, in
general, there is no a priori small parameter that justifies such a truncation
and the validity of the truncation must be considered on a case-by-case basis.
We will return to this issue in Sec. IV.5.

By repeating these steps for the von Neumann equation, we obtain a Fokker-
Planck equation for the FWF $\chi({\bm{\theta}},{\bm{\phi}},t)$, resulting for
the case of the Ising Hamiltonian in the set of ordinary differential
equations for the spin angles

$\displaystyle\frac{d}{dt}\theta_{n}(t)=$ $\displaystyle 0,$ (43a)
$\displaystyle\frac{d}{dt}\phi_{n}(t)=$ $\displaystyle\sqrt{3}\sum_{m\neq
n}J_{mn}\cos\theta_{m}(t).$ (43b)

Transforming to Cartesian coordinates yields

$\displaystyle\frac{d}{dt}\begin{pmatrix}s_{n}^{x}\\\ s_{n}^{y}\\\
s_{n}^{z}\end{pmatrix}=\sum_{m\neq
n}J_{mn}s_{m}^{z}\begin{pmatrix}-s_{n}^{y}\\\ +s_{n}^{x}\\\ 0\end{pmatrix},$
(44a)

which are precisely the mean-field EOMs for the Ising Hamiltonian in DTWA.

###  IV.3 Dephasing, decay, and incoherent pump

The coupling of spins to Markovian reservoirs can be straightforwardly
included in our hybrid approach. Applying the mappings between the Hilbert
space and phase space, Eq. (33), we have for the dephasing Lindbladian

$\displaystyle\hat{\sigma}^{z}\hat{\rho}\hat{\sigma}^{z}-\hat{\rho}\enspace\longleftrightarrow\enspace
2\frac{\partial^{2}}{\partial\phi^{2}}\chi.$ (45)

Using a change of variables, we can transform the Fokker-Planck equation
resulting from Eq. (45) to Cartesian coordinates, obtaining the set of SDEs

$\displaystyle ds_{x}=$ $\displaystyle-2s_{x}dt-2s_{y}dW_{\phi},$ (46a)
$\displaystyle ds_{y}=$ $\displaystyle-2s_{y}dt+2s_{x}dW_{\phi},$ (46b)
$\displaystyle ds_{z}=$ $\displaystyle\,0.$ (46c)

These equations reproduce Eqs. (13)-(15) of Ref. [22] if we make the
substitution $t\rightarrow 2t/\Gamma_{\phi}$.

In a similar way, we can derive the mappings for incoherent gains and losses,

$\displaystyle\hat{\sigma}^{\pm}\hat{\rho}\hat{\sigma}^{\mp}-\frac{1}{2}\\{\hat{\sigma}^{\mp}\hat{\sigma}^{\pm},\hat{\rho}\\}\,\longleftrightarrow\,-\frac{\partial}{\partial\theta}\left(\cot\theta\pm\frac{\csc\theta}{\sqrt{3}}\right)\chi$
$\displaystyle+\frac{1}{2}\frac{\partial^{2}}{\partial\phi^{2}}\left(1+2\cot^{2}\theta\pm\frac{2\cot\theta\csc\theta}{\sqrt{3}}\right)\chi.$
(47)

Changing to Cartesian coordinates similarly reproduces the deterministic parts
given by Eqs. (24) of Ref. [22]. One finds, however, that the diffusion matrix
in Cartesian coordinates is not positive semidefinite and no corresponding SDE
exists. In stark contrast to this, the parametrization with respect to
$\theta,\phi$ does have a corresponding set of SDEs that can immediately be
deduced from Eq. (47) as

$\displaystyle d\theta=$
$\displaystyle\left(\cot\theta\pm\frac{\csc\theta}{\sqrt{3}}\right)dt,$ (48a)
$\displaystyle d\phi=$
$\displaystyle\sqrt{1+2\cot^{2}\theta\pm\frac{2\cot\theta\csc\theta}{\sqrt{3}}}\,\,dW.$
(48b)

Thus the hybrid discrete-continuous approach can incorporate incoherent
processes which lead to classical noise terms in the EOMs.

###  IV.4 Gauge freedom of SU(2) Wigner functions and sampling of initial
states

We have seen in Sec. III that all pure spin states have a continuous SU(2)
Wigner function $W(\theta,\phi)$ which is non-positive and thus cannot be
sampled by Monte Carlo methods. This complicates explicit averaging over the
initial Wigner distribution. This problem can be overcome by relating the
continuous and discrete Wigner representations of spins using a gauge freedom
[35].

The SU(2) phase-point operators of Eq. (21) are orthogonal to all functions
that only contain spherical harmonics $\text{Y}_{l,m}$ with $l\geq 2$,

$f(\theta,\phi)=\sum_{l=2}^{\infty}\sum_{m=-l}^{l}C_{l,m}\text{Y}_{l,m}(\theta,\phi),$
(49)

which means that

$\int d\Omega\,\hat{A}(\theta,\phi)\,f(\theta,\phi)=0.$ (50)

The function $f$ can be chosen real by setting $C_{l,-m}=(-1)^{m}C_{l,m}$.
Thus the representation of a spin state $\hat{\rho}$ in terms of a SU(2)
Wigner function is not unique but has a gauge freedom

$W(\theta,\phi)\,\equiv\,W(\theta,\phi)+f(\theta,\phi).$ (51)

The initial discrete Wigner coefficients $W_{\bm{\alpha}0}$ are positive in
many important cases. We now argue that the initial FWF
$\chi_{0}(\theta,\phi)$ can be expressed in terms of $W_{\bm{\alpha}0}$ as

$\displaystyle\chi_{0}(\theta,\phi)=$
$\displaystyle\sum_{\bm{\alpha}}\delta(\theta-\theta_{\bm{\alpha}})\delta(\phi-\phi_{\bm{\alpha}})W_{\bm{\alpha}0},$
(52)

where the initial phase point operators
$\hat{A}_{\bm{\alpha}}=\hat{A}(\theta_{\bm{\alpha}},\phi_{\bm{\alpha}})$
correspond to angles

$\displaystyle\theta_{00}$
$\displaystyle=\theta_{01}=\pi-\arccos{\frac{1}{\sqrt{3}}},$
$\displaystyle\theta_{10}$
$\displaystyle=\theta_{11}=\arccos{\frac{1}{\sqrt{3}}},$ (53)
$\displaystyle\phi_{00}$
$\displaystyle=\frac{7\pi}{4},\,\phi_{01}=\frac{3\pi}{4},\,\phi_{10}=\frac{\pi}{4},\,\phi_{11}=\frac{5\pi}{4}.$

Substituting this into Eq. (25a) reproduces the initial density operator. As
an example, the spin down state is sampled from two points as

$\displaystyle\chi^{(2p)}(\theta,\phi)=$
$\displaystyle\delta\left(\theta-\arccos\frac{1}{\sqrt{3}}\right)$
$\displaystyle\cdot\frac{1}{2}\left[\delta\left(\phi-\frac{\pi}{4}\right)+\delta\left(\phi-\frac{5\pi}{4}\right)\right].$
(54)

Next, it follows from the completeness relation of spherical harmonics,

$\displaystyle\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\text{Y}^{*}_{lm}(\theta_{\bm{\alpha}},\phi_{\bm{\alpha}})\text{Y}_{lm}(\theta,\phi)$
$\displaystyle\qquad=\delta(\phi-\phi_{\bm{\alpha}})\,\delta(\cos\theta-\cos\theta_{\bm{\alpha}}),$

that the difference between $\chi_{0}(\theta,\phi)$ of Eq. (52) and a directly
evaluated FWF pertaining to the initial state contains only components of type
$\sin(\theta)f(\theta,\phi)/2\pi$.

Using Eq. (52) and the fact that all eigenstates of the Pauli matrices have
positive discrete Wigner coefficients $W_{\bm{\alpha}0}$, we can sample the
initial phase space distribution from a discrete set of points. Using the
discrete set of points decreases furthermore the statistical error when a
finite number of samples is taken. Since all pure states (spin coherent
states) are connected by unitary operations, we can generate positive discrete
distributions for them as well. Finally, every mixed state can be represented
as a mixture of pure states and therefore every single particle state is
accessible by a classical Monte Carlo sampling.

###  IV.5 Validity of the interaction truncation in DTWA

The original formulation of the DTWA [10] does not contain information on the
range of validity of the interaction truncation. Instead, the quality of this
approximation has only been characterized empirically, by comparing the
predictions of the DTWA with exact (numerical) solutions for a given
Hamiltonian. It was found in Ref. [10] that the DTWA results reproduce the
short time dynamics of macroscopic spin observables well if the interaction
couples many spins in a similar way. We now provide an explanation for this
empirical observation by deriving conditions under which the truncation
approximation is justified.

Our aim is to show that the second and third derivatives in Eq. (42) generated
by the two-body interactions of Eq. (41) can indeed be neglected in the case
of macroscopic collective spin dynamics or for short interaction times. For
simplicity, we assume the extreme limit of $J_{mn}=J$, i.e. consider all-to-
all interactions. We then find that the FWF satisfies a PDE of the form

$\displaystyle\frac{\partial}{\partial t}\chi=$
$\displaystyle-\sum_{n}\frac{\partial}{\partial\phi_{n}}(A_{n}\chi)+\sum_{n}\sum_{m\neq
n}\frac{\partial^{2}}{\partial\theta_{m}\partial\phi_{n}}(D_{mn}\chi)$
$\displaystyle-\sum_{n}\sum_{m\neq
n}\frac{\partial^{3}}{\partial\phi_{m}^{2}\partial\phi_{n}}(G_{mn}\chi),$ (55)

where the coefficients are given by

$\displaystyle A_{n}=$ $\displaystyle\sqrt{3}J\sum_{m\neq n}\cos\theta_{m},$
(56a) $\displaystyle D_{mn}=$ $\displaystyle
J\frac{2\csc\theta_{m}-3\sin\theta_{m}}{\sqrt{3}},$ (56b) $\displaystyle
G_{mn}=$ $\displaystyle 2J\frac{\cot\theta_{m}\csc\theta_{m}}{\sqrt{3}}.$
(56c)

If all spins are aligned, $\theta_{i}\approx\theta_{j}\;\forall\;i,j$, then
all terms in the sum of Eq. (56a) interfere constructively. The drift
coefficients are therefore of order $A_{n}\sim O(N)$, while the other
coefficients are of order $D_{mn},G_{mn}\sim O(1)$. By introducing scaled
$\phi$-variables $\tilde{\phi}_{n}=\phi_{n}/N$, we obtain new coefficients

$\displaystyle\tilde{A}_{n}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}J}{N}\sum_{m\neq n}\cos\theta_{m}\,$
$\displaystyle=\,O(1)$ $\displaystyle\tilde{D}_{mn}$ $\displaystyle=$
$\displaystyle\frac{J}{N}\frac{2\csc\theta_{m}-3\sin\theta_{m}}{\sqrt{3}}\,$
$\displaystyle=\,O(N^{-1}),$ (57) $\displaystyle\tilde{G}_{mn}$
$\displaystyle=$
$\displaystyle\frac{J}{N^{3}}\frac{2\cot\theta_{m}\csc\theta_{m}}{\sqrt{3}}\,$
$\displaystyle=\,O(N^{-3}).$

The sums over all first derivatives and second derivatives each scale as
$O(N)$, but the sum over the third derivatives scales as $O(N^{-1})$. Hence,
for large values of $N\gg 1$ we can neglect the third derivatives and obtain a
generalized Fokker-Planck equation. Note however that the diffusion matrix
associated with mixed derivatives is not positive definite. To arrive at the
deterministic equations of the DTWA, all derivatives higher than first-order
must be disregarded. To deduce the conditions under which this is justified,
we perform a change of variables using the mean
$\theta=\frac{1}{N}\sum_{i}\theta_{i}$ and difference
$\delta_{n}=\theta_{n}-\theta_{n+1}$ angles which yields a new generalized FPE

$\displaystyle\frac{\partial}{\partial
t}\tilde{\chi}=-\sum_{n}\frac{\partial}{\partial\tilde{\phi}_{n}}\left(\tilde{A}_{n}\tilde{\chi}\right)+\sum_{n}\frac{\partial^{2}}{\partial\tilde{\phi}_{n}\partial\theta}\left(\frac{1}{N}\tilde{D}_{mn}\tilde{\chi}\right)$
$\displaystyle+\sum_{mn}\frac{\partial^{2}}{\partial\tilde{\phi}_{m}\partial\delta_{n}}\left((\tilde{D}_{mn+1}-\tilde{D}_{mn})\tilde{\chi}\right).$
(58)

The second term on the r.h.s. is $O(1)$ and can therefore be neglected in
comparison to the first term which is $O(N)$. The third term is, however, also
of order $O(N)$, but contains only differences of polar angles

$\displaystyle\tilde{D}_{mn+1}-\tilde{D}_{mn}=$ (59)
$\displaystyle\quad\frac{J}{N}\frac{3(\sin\theta_{n+1}-\sin\theta_{n})-2(\csc\theta_{n+1}-\csc\theta_{n})}{\sqrt{3}}.$

If, as assumed, there is an all-to-all spin coupling and initially all the
spins are aligned with high probability,
$\theta_{n+1}(0)\approx\theta_{n}(0)$, then the difference of the diffusion
coefficients is small and the corresponding term in Eq. (58) can be neglected.

This explains the findings in Ref. [10], where collective observables are well
reproduced by the DTWA if the interaction Hamiltonian has a high effective
coordination number. One can similarly explain why the DTWA generally gives
good predictions for short-time dynamics: At small times, the dynamical
variables $\theta(t)$ and $\phi(t)$ have values of Eq. (53) with probabilities
given by the discrete Wigner coefficients $W_{\bm{\alpha}}$ of the initial
state. For these values the diffusion coefficients $D_{mn}$ vanish
identically.

##  V Dynamics of a laser-driven array of Rydberg atoms

We illustrate the performance of the DCTWA by applying it to an experimentally
relevant system of driven, dissipative, interacting spins. Specifically, we
consider an array of atoms driven by a resonant laser to the strongly
interacting Rydberg state [36]. Denoting the atomic ground state by
$|\negmedspace\downarrow\rangle$ and the excited Rydberg state by
$|\negmedspace\uparrow\rangle$, the Hamiltonian of the system is given by

$\displaystyle\hat{H}=\Omega\sum_{n}\hat{\sigma}_{n}^{x}+\frac{1}{2}\sum_{m\neq
n}\frac{J}{|m-n|^{\alpha}}\hat{\sigma}_{m}^{rr}\hat{\sigma}_{n}^{rr},$ (60)

where $\Omega$ is the Rabi frequency of the resonant laser,
$\hat{\sigma}_{n}^{rr}=(\hat{\mathbbm{1}}+\hat{\sigma}_{n}^{z})/2$ is the
projector onto the Rydberg state, $J$ is the interaction strength and $\alpha$
determines the interaction range (e.g., $\alpha=6$ for van der Waals
interactions). We include a local dephasing with rate $\kappa$ and incoherent
decay (deexcitation) of the Rydberg state with rate $\gamma$ via

$\displaystyle\hat{L}_{n}^{\kappa}=\sqrt{\kappa}\hat{\sigma}_{n}^{z},$
$\displaystyle\quad\hat{L}_{n}^{\gamma}=\sqrt{\gamma}\hat{\sigma}_{n}^{-},$
(61)

where
$\hat{\sigma}_{n}^{\pm}=i(\hat{\sigma}_{n}^{x}\pm\hat{\sigma}_{n}^{y})/2$.

Without interactions ($J=0$), the laser field induces damped Rabi oscillations
of all the atoms, resulting in a stationary state with

$\langle\hat{S}_{z}\rangle=\frac{1}{N}\sum_{n}\langle\hat{\sigma}^{z}_{n}\rangle\to\langle\hat{\sigma}^{z}\rangle=-\frac{\gamma\kappa+(\gamma/2)^{2}}{2\Omega^{2}+\gamma\kappa+(\gamma/2)^{2}}.$

In the presence of interactions, $J\neq 0$, an atom in the Rydberg state
shifts the Rydberg transition of all the surrounding atoms out of resonance.
Within the blockade distance, this shift is sufficiently large to suppress the
excitation of the other atoms [36]. As a consequence, the steady-state value
of the effective spin polarization $\langle\hat{S}_{z}\rangle$ is reduced.

Mapping the Lindbladian to the spin phase space and truncating interaction
contributions leads to the set of SDEs

$\displaystyle d\theta_{n}=$
$\displaystyle\left[-2\Omega\sin\phi_{n}+\gamma\left(\cot\theta_{n}-\frac{\csc\theta_{n}}{\sqrt{3}}\right)\right]dt,$
(62a) $\displaystyle d\phi_{n}=$
$\displaystyle-\bigg{(}2\Omega\cot\theta_{n}\cos\phi_{n}+\frac{J}{2}\sum_{m\neq
n}\frac{1-\sqrt{3}\cos\theta_{m}}{|m-n|^{\alpha}}\bigg{)}dt$
$\displaystyle+\sqrt{\gamma\left(1+2\cot^{2}\theta_{n}-\frac{2\cot\theta_{n}\csc\theta_{n}}{\sqrt{3}}\right)+4\kappa}\,dW_{\phi_{n}}.$
(62b)

We evaluated these equations for atoms in a one-dimensional lattice with
periodic boundary conditions using the DifferentialEquations package [37] of
the Julia programming language [38]. The observables were calculated by
averaging $92\times 10^{3}$ trajectories, and we have verified the convergence
of the calculations by further increasing the number of trajectories and
observing no significant variation of the results. We note that the first
moments only require $\sim 10^{3}$ trajectories to convergence properly. Due
to the dissipative fluctuations, correlations and other higher moments in the
steady state require many more trajectories.

Figure 2: First collective moments
$\hat{S}_{\mu}=\frac{1}{N}\sum_{n=1}^{N}\hat{\sigma}_{n}^{\mu}$ for $N=10$
atoms in a one-dimensional lattice (periodic boundary conditions) excited to
the Rydberg state by a resonant laser, as obtained from DCTWA (solid blue
lines), exact solution of the master equation (dashed-dotted black line), and
the mean-field calculations (solid orange lines in the insets). The parameters
are $\Omega=0.3J$, $\alpha=6$ and $\gamma=\kappa=10^{-2}J$. Figure 3: First
collective moments
$\hat{S}_{\mu}=\frac{1}{N}\sum_{n=1}^{N}\hat{\sigma}_{n}^{\mu}$ with
parameters as given in Fig. 2, but with initial state sampling according to
Eqs. (54), (63), and (64) (solid lines) and exact results (dashed-dotted black
line). The results obtained from the 4p-sampling (orange solid line) and
$\infty$p-sampling (green solid line) coincide.

In Fig. 2 we show the time dependence of the average spin polarizations
$\langle\hat{S}_{\mu}(t)\rangle=\frac{1}{N}\sum_{n=1}^{N}\langle\hat{\sigma}_{n}^{\mu}(t)\rangle$
along the $x,y,z$ directions for $N=10$ atoms subject to relatively large Rabi
frequency $\Omega=0.3J$ and weak damping and dephasing rates
$\kappa=\gamma=0.01J$, as obtained via our hybrid DCTWA (blue solid lines).
For comparison, we also show the results obtained from exact solutions of the
density matrix equations under the same conditions (black dashed-dotted
lines). We observe that the results of DCTWA are in good agreement with the
exact simulations, especially for the final stationary values of the spin
polarizations. In the insets of Fig. 2 we also show the results of the mean-
field calculations (orange solid lines). The mean-field equations do not
contain the stochastic terms resulting from decay and dephasing, and no
averaging over an initial distribution is performed. As a consequence, the
persistence of Rabi oscillations is substantially overestimated and also the
stationary values of the collective spin deviate from the exact results.

In Eqs. (7) we have introduced a discrete representation of the spin down
state which uses two of the four discrete phase space points (2$p$ sampling).
According to Eq. (5) this means that the signs of the Cartesian $x$\- and
$y$-component are strictly anti-correlated. While this correctly reproduces
the initial density operator, time evolving these anti-correlations according
to the _approximate_ dynamics of the DTWA can affect the dynamics of
interacting spins [39, 17]. Let us, therefore, employ a sampling scheme of
initial states that uses more discrete phase points as this reduces the
relevance of such correlations, while still faithfully representing the
initial state. Specifically, consider the set of four discrete phase points
that is generated by applying a $\pi/2$ rotation around the $z$-axis. By
combining both sets of four points each, we can express any given state using
eight points. The spin down state, e.g., is described by four points with non-
vanishing Wigner coefficients

$\displaystyle\chi^{(4p)}(\theta,\phi)=$
$\displaystyle\delta\left(\theta-\arccos\frac{1}{\sqrt{3}}\right)$
$\displaystyle\cdot\frac{1}{4}\sum_{n=1}^{4}\delta\left(\phi-\frac{(2n-1)\pi}{4}\right),$
(63)

which corresponds to setting the Cartesian component $s_{z}=-1$ and
independently drawing $\pm 1$ with equal probability for the $x$ and $y$
component. Similarly we can consider a sampling from an even larger number of
discrete phase points arriving eventually at a continuous distribution arising
from rotations around the $z$ axis. This results in

$\displaystyle\chi^{(\infty p)}(\theta,\phi)=$
$\displaystyle\frac{1}{2\pi}\delta\left(\theta-\arccos\frac{1}{\sqrt{3}}\right).$
(64)

The dynamics produced from these sampling schemes are compared in Fig. 3. The
2$p$ sampling corresponds to the blue lines of Fig. 2 and differs from the
4$p$\- and $\infty p$ samplings which produce identical dynamics that resemble
the exact results much more closely. The drastic initial overshooting in
$\langle\hat{S}_{x}\rangle$ of the 2p-method is corrected. Furthermore, the
extremal points of all components are slightly shifted in time which
significantly increases the quantitative agreement with the exact data at
short times as well as the qualitative agreement at short-to-intermediate time
scales. The improved agreement with exact results can be understood as
follows: The DTWA amounts to neglecting the cross diffusion terms in the
second line of Eq. (58), which do have the same scaling ${\cal O}(N)$ with the
number $N$ of spins as the first term in Eq. (58) and thus are not a priori
small. In the $\infty$p-sampling scheme discussed above the initial Wigner
distribution is however homogeneous in all $\tilde{\phi}_{n}$ and thus the
derivative
$(\partial^{2}/\partial\tilde{\phi}_{n}\partial\delta_{n})\tilde{\chi}(t=0)$
vanishes.

Figure 4: Stationary Rydberg-Rydberg correlations of atoms, as obtained from
the DCTWA and exact solution of the density matrix equations, for the same
parameters as in Fig. 2. The correlations are evaluated at $Jt=200$ when the
steady state has long been reached.

In Fig. 4 we finally show the steady-state correlations of Rydberg excitations
of the atoms in the lattice. The competition between the resonant laser
excitation and Rydberg blockade of nearest-neighbor atoms leads to a density-
wave of Rydberg excitations. In one dimension, the resulting steady-state
correlations always decay exponentially, while the correlation length is very
small for the two-level driving scheme considered here [40]. We note that the
correlations are well captured by our DCTWA method, while they cannot be
calculated by using a simple mean-field approach.

##  VI Summary

We have presented a practical approach for a semiclassical description of
interacting spins which is a hybrid discrete-continuous generalization of the
discrete truncated Wigner approximation. In our DCTWA, the quantum state of
individual spins is represented by a Wigner function in a continuous rather
than a discrete phase space, and, as in standard DTWA, interactions are
treated via mean-field factorization. Quantum fluctuations are taken into
account in lowest order by averaging over the quasiprobability-distribution of
initial states. The advantage of a Wigner representation in a continuous phase
space is that the corresponding equation of motion is a partial differential
equation. Under specific, quantifiable conditions, this equation can be
approximated by a Fokker-Planck equation with positive definite diffusion and
mean-field interaction contributions, which can then efficiently be simulated
by solving ordinary SDEs for the spins using an angular representation. An
important property of the continuous representation of spin states is the
overcompleteness of the corresponding phase-point operators which leads to a
gauge freedom in the continuous Wigner function. This gauge freedom can be
used to overcome the main drawback of a continuous Wigner function, namely its
non-positivity for pure spin states. Exploiting this freedom and mapping the
continuous Wigner function of typical initial states to their discrete
counterpart allows an averaging over the initial state by Monte Carlo
sampling.

The DCTWA allows for a rigorous derivation of the truncation approximation and
yields conditions for its applicability. Hence, we were able to explain the
empirically observed range of validity of the DTWA. Furthermore, the DCTWA
allows us to include Markovian reservoir couplings leading to dephasing, decay
or incoherent pumping in a straightforward way. Disregarding the noise terms
in the SDEs, resulting from reservoir couplings and transforming to Cartesian
spin coordinates, we reproduce the standard DTWA equations [10]. Considering
dephasing reproduces the stochastic equations of Ref. [22]. Decay and
incoherent pumping, on the other hand, lead to non-classical noise terms in
the equations of motion of Cartesian spin components used in the standard
DTWA. These processes can only be treated by stochastic simulations in an
angular representation in which the diffusion terms are positive definite.

We have illustrated the performance of the DCTWA by considering a small one-
dimensional array of atoms, resonantly driven into a Rydberg state in the weak
damping regime under conditions of a nearest-neighbor Rydberg blockade.
Comparison of the time dependence of collective spin observables and steady-
state spin-spin correlations showed very good agreement with exact
simulations, in contrast to the mean-field calculations.

Our approach paves the way for systematic improvements of the standard DTWA,
which will be the subject of future work.

###### Acknowledgements.

We would like to thank A.-M. Rey, J. Schachenmayer, and H. Weimer for fruitful
discussions. We thank J. Schachenmayer for pointing the benefit of an initial-
state sampling scheme with a larger number of phase points out to us. We
acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) via the Collaborative Research Center SFB/TR185
(Project No. 277625399) and the priority programm SPP 1929 (Project No.
273920612) and of the Alexander von Humboldt Foundation via the Research Group
Linkage Programme. D.P. was also supported by the EU QuantERA Project PACE-IN
(GSRT Grant No. T11EPA4-00015).

##  Appendix A CLASSICAL SPIN EQUATIONS OF MOTION

The proof of Eq. (10) follows from projecting $\frac{d}{dt}\hat{A}(t)$ of Eq.
(8) onto a single element $s_{j}(t)$

$\displaystyle\dot{s}_{j}=$
$\displaystyle\text{Tr}\\{\hat{\sigma}_{j}\frac{d}{dt}\hat{A}\\}$
$\displaystyle=$
$\displaystyle-\frac{i}{2}\sum_{kl}s_{l}h_{k}\text{Tr}\\{\hat{\sigma}_{j}[\hat{\sigma}_{k},\hat{\sigma}_{l}]\\}$
$\displaystyle=$ $\displaystyle
2\sum_{kl}\epsilon_{klm}s_{l}h_{k}\text{Tr}\\{\hat{\sigma}_{j}\hat{\sigma}_{m}\\}$
$\displaystyle=$ $\displaystyle
2\sum_{kl}\epsilon_{jkl}s_{l}\frac{\partial{\cal H}}{\partial s_{k}}.$ (65)

In the last line we have used the fact that the classical Hamilton function
${\cal H}=\sum_{k}h_{k}s_{k}$ is always linear in the spin $\bm{s}$ and
therefore $h_{k}=\frac{\partial{\cal H}}{\partial s_{k}}$.

##  Appendix B PHASE SPACE MAPPINGS

As demonstrated in Sec. IV.1, an operator acting on a state $\hat{\rho}$ can
be translated into a differential operator acting on the flattened Wigner
function $\chi(\theta,\phi,t)=\sin(\theta)W(\theta,\phi)/2\pi$. Repeating the
steps performed for $\hat{\sigma}^{z}\hat{\rho}$ for all Pauli matrices acting
from the left and right yields

$\displaystyle\hat{\sigma}^{x}\hat{\rho}\leftrightarrow$
$\displaystyle\left[\sqrt{3}\sin\theta\cos\phi-\frac{\partial}{\partial\theta}\left(\sqrt{3}\cos\theta\cos\phi-i\sin\phi\right)+\frac{\partial}{\partial\phi}\left(\frac{\csc\theta\sin\phi}{\sqrt{3}}+i\cot\theta\cos\phi\right)+\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\csc\theta\cos\phi}{\sqrt{3}}\right]\chi,$
(66a) $\displaystyle\hat{\rho}\hat{\sigma}^{x}\leftrightarrow$
$\displaystyle\left[\sqrt{3}\sin\theta\cos\phi-\frac{\partial}{\partial\theta}\left(\sqrt{3}\cos\theta\cos\phi+i\sin\phi\right)+\frac{\partial}{\partial\phi}\left(\frac{\csc\theta\sin\phi}{\sqrt{3}}-i\cot\theta\cos\phi\right)+\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\csc\theta\cos\phi}{\sqrt{3}}\right]\chi,$
(66b) $\displaystyle\hat{\sigma}^{y}\hat{\rho}\leftrightarrow$
$\displaystyle\left[-\sqrt{3}\sin\theta\sin\phi+\frac{\partial}{\partial\theta}\left(\sqrt{3}\cos\theta\sin\phi+i\cos\phi\right)+\frac{\partial}{\partial\phi}\left(\frac{\csc\theta\cos\phi}{\sqrt{3}}-i\cot\theta\sin\phi\right)-\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\csc\theta\sin\phi}{\sqrt{3}}\right]\chi,$
(66c) $\displaystyle\hat{\rho}\hat{\sigma}^{y}\leftrightarrow$
$\displaystyle\left[-\sqrt{3}\sin\theta\sin\phi+\frac{\partial}{\partial\theta}\left(\sqrt{3}\cos\theta\sin\phi-i\cos\phi\right)+\frac{\partial}{\partial\phi}\left(\frac{\csc\theta\cos\phi}{\sqrt{3}}+i\cot\theta\sin\phi\right)-\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\csc\theta\sin\phi}{\sqrt{3}}\right]\chi,$
(66d) $\displaystyle\hat{\sigma}^{z}\hat{\rho}\leftrightarrow$
$\displaystyle\left[-\sqrt{3}\cos\theta-\frac{\partial}{\partial\theta}\frac{3\sin\theta-2\csc\theta}{\sqrt{3}}+\frac{\partial}{\partial\phi}i-\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\cot\theta\csc\theta}{\sqrt{3}}\right]\chi,$
(66e) $\displaystyle\hat{\rho}\hat{\sigma}^{z}\leftrightarrow$
$\displaystyle\left[-\sqrt{3}\cos\theta-\frac{\partial}{\partial\theta}\frac{3\sin\theta-2\csc\theta}{\sqrt{3}}-\frac{\partial}{\partial\phi}i-\frac{\partial^{2}}{\partial\phi^{2}}\frac{2\cot\theta\csc\theta}{\sqrt{3}}\right]\chi.$
(66f)

Any differential equation with respect to $\hat{\rho}$ in Hilbert space can
thus be translated into a PDE in the phase space.

##  Appendix C COMPLEX STEREOGRAPHIC PROJECTION MAPPINGS

The parametrization of the phase-point operator $\hat{A}(\Omega)$ is not
uniquely given by a pair of angles $\theta,\phi$. Instead, one can introduce a
steoreographic projection onto a complex plane with coordinates
$\beta\in\mathbb{C}$. In this case, Eq. (22) can be equivalently expressed as

$\displaystyle\bm{s}(\beta)=\frac{\sqrt{3}}{1+|\beta|^{2}}\Bigl{(}\beta+\beta^{*},-i(\beta-\beta^{*}),-1+|\beta|^{2}\Bigl{)}^{T}.$
(67)

The formulation with respect to the angles $\theta,\phi$ is recovered by
substituting $\beta=\tan(\theta/2)e^{-i\phi}$. The new integral measure is
given by

$\displaystyle\int d\Omega=\int d^{2}\beta\,\frac{2}{\pi(1+|\beta|^{2})^{2}},$
(68)

such that $1=\int
d^{2}\beta\,\frac{2}{\pi(1+|\beta|^{2})^{2}}W(\beta^{*},\beta)$. The matrices
$\hat{A},\frac{\partial\hat{A}}{\partial\beta},\frac{\partial\hat{A}}{\partial\beta^{*}},\frac{\partial^{2}\hat{A}}{\partial\beta^{*}\partial\beta}$
span the Hilbert space, where $\beta$ and $\beta^{*}$ are treated as
independent variables. The discrete phase-point operators
$\hat{A}_{\bm{\alpha}}=\hat{A}(\beta_{\bm{\alpha}}^{*},\beta_{\bm{\alpha}})$
are given at points

$\displaystyle\beta_{\bm{\alpha}}=(-1)^{\alpha_{1}+\alpha_{2}}(-1)^{\frac{1+2\alpha_{1}}{4}}\frac{(-1)^{\alpha_{1}}+\sqrt{3}}{\sqrt{2}}.$
(69)

We can similarly derive a set of mappings for the new FWF
$\chi(\beta^{*},\beta)=\frac{2}{\pi(1+|\beta|^{2})^{2}}W(\beta^{*},\beta)$

$\displaystyle\hat{\sigma}^{x}\hat{\rho}\leftrightarrow$
$\displaystyle\left[\sqrt{3}\frac{\beta+\beta^{*}}{1+|\beta|^{2}}-\frac{\partial}{\partial\beta}\frac{c_{+}}{6}(1-\beta^{2})-\frac{\partial}{\partial\beta^{*}}\frac{c_{-}}{6}(1-\beta^{*2})+\frac{\partial^{2}}{\partial\beta^{*}\partial\beta}\frac{(\beta+\beta^{*})(1+|\beta|^{2})}{\sqrt{3}}\right]\chi(\beta^{*},\beta),$
(70a) $\displaystyle\hat{\rho}\hat{\sigma}^{x}\leftrightarrow$
$\displaystyle\left[\sqrt{3}\frac{\beta+\beta^{*}}{1+|\beta|^{2}}-\frac{\partial}{\partial\beta}\frac{c_{-}}{6}(1-\beta^{2})-\frac{\partial}{\partial\beta^{*}}\frac{c_{+}}{6}(1-\beta^{*2})+\frac{\partial^{2}}{\partial\beta^{*}\partial\beta}\frac{(\beta+\beta^{*})(1+|\beta|^{2})}{\sqrt{3}}\right]\chi(\beta^{*},\beta),$
(70b) $\displaystyle\hat{\sigma}^{y}\hat{\rho}\leftrightarrow$
$\displaystyle\left[-i\sqrt{3}\frac{\beta-\beta^{*}}{1+|\beta|^{2}}-\frac{\partial}{\partial\beta}\frac{ic_{+}}{6}(1+\beta^{2})+\frac{\partial}{\partial\beta^{*}}\frac{ic_{-}}{6}(1+\beta^{*2})-\frac{\partial^{2}}{\partial\beta^{*}\partial\beta}i\frac{(\beta-\beta^{*})(1+|\beta|^{2})}{\sqrt{3}}\right]\chi(\beta^{*},\beta),$
(70c) $\displaystyle\hat{\rho}\hat{\sigma}^{y}\leftrightarrow$
$\displaystyle\left[-i\sqrt{3}\frac{\beta-\beta^{*}}{1+|\beta|^{2}}-\frac{\partial}{\partial\beta}\frac{ic_{-}}{6}(1+\beta^{2})+\frac{\partial}{\partial\beta^{*}}\frac{ic_{+}}{6}(1+\beta^{*2})-\frac{\partial^{2}}{\partial\beta^{*}\partial\beta}i\frac{(\beta-\beta^{*})(1+|\beta|^{2})}{\sqrt{3}}\right]\chi(\beta^{*},\beta),$
(70d) $\displaystyle\hat{\sigma}^{z}\hat{\rho}\leftrightarrow$
$\displaystyle\left[\sqrt{3}\frac{1-|\beta|^{2}}{1+|\beta|^{2}}+\frac{\partial}{\partial\beta}\frac{c_{+}}{3}\beta+\frac{\partial}{\partial\beta^{*}}\frac{c_{-}}{3}\beta^{*}+\frac{\partial^{2}}{\partial\beta^{*}\partial\beta}\frac{1-|\beta|^{4}}{\sqrt{3}}\right]\chi(\beta^{*},\beta),$
(70e) $\displaystyle\hat{\rho}\hat{\sigma}^{z}\leftrightarrow$
$\displaystyle\left[\sqrt{3}\frac{1-|\beta|^{2}}{1+|\beta|^{2}}+\frac{\partial}{\partial\beta}\frac{c_{-}}{3}\beta+\frac{\partial}{\partial\beta^{*}}\frac{c_{+}}{3}\beta^{*}+\frac{\partial^{2}}{\partial\beta^{*}\partial\beta}\frac{1-|\beta|^{4}}{\sqrt{3}}\right]\chi(\beta^{*},\beta),$
(70f)

where $c_{\pm}=\sqrt{3}\pm 3$.

##  Appendix D BENCHMARKING OF THE DCTWA

We show that the DCTWA precisely reproduces the exact dynamics of a single
spin, driven by an external field and subject to spontaneous decay. In Fig. 5
we plot the results of the DCTWA simulation of
$\langle\hat{\sigma}^{z}(t)\rangle$ and compare them to exact results for a
spin initially in the eigenstate of $\hat{\sigma}^{y}$, $\hat{\sigma}^{z}$
with eigenvalue $-1$ and evolving according to the Hamiltonian
$\hat{H}=\frac{g}{2}\hat{\sigma}^{x}$ and Lindblad generator
$\hat{L}=\sqrt{\gamma}\hat{\sigma}^{-}$. We choose parameters identical to
that in Ref. [24] and find perfect agreement with the exact results, both
dynamically as well as in the long time limit $\gamma t=15$.

Figure 5: Time evolution of $\langle\hat{\sigma}^{z}\rangle$ with initial
states being polarized along the negative Z/Y axis obtained from numerical
integration of the DCTWA (blue/orange solid lines respectively) and exact
results (black lines). The inset shows the long-time limit (orange crosses)
for varying $g/\gamma$ and exact steady state results (black line).

A similar excellent agreement was obtained using the OSDTWA (open system
discrete truncated Wigner approximation) introduced in Ref. [24], which leads
to a simpler, but larger set of four dynamical equations for every spin,
describing its Cartesian components and one additional degree of freedom
$S^{0}$ describing the norm:

$\displaystyle\dot{S}^{x}=$ $\displaystyle-\frac{\gamma}{2}S^{x},$ (71)
$\displaystyle\dot{S}^{y}=$ $\displaystyle-\frac{\gamma}{2}S^{y}-gS^{z},$ (72)
$\displaystyle\dot{S}^{z}=$
$\displaystyle-\frac{\gamma}{2}(S^{0}+S^{z})+gS^{y},$ (73)
$\displaystyle\dot{S}^{0}=$ $\displaystyle-\frac{\gamma}{2}(S^{0}+S^{z})$ (74)

In the OSDTWA the spin dynamics are obtained by evolving these equations from
initial values corresponding to the discrete Wigner distribution of the
initial state as well as $S^{0}(0)=1$ and probabilistically performing quantum
jumps. Whether or not a jump occurs at a given time step is determined by the
jump probability $\delta p=\frac{\gamma}{2}(S^{0}+S^{z})$. We note that
although this approach leads to a perfect agreement with exact simulations for
a single spin initially prepared in the state $|\downarrow\rangle$, it may
lead to some artifacts in other cases that must be corrected by hand. Let us
consider, e.g., as an initial condition the eigenstate of $\hat{\sigma}^{y}$
with eigenvalue $-1$, which corresponds to realizations starting with
$S^{y}(0)=-1$, $S^{0}(0)=1$ and $S^{x}(0)$, $S^{z}(0)$ having an equal
probability of being $\pm 1$. More specifically we choose the trajectory with
$S^{z}(0)=-1$. Initially, the quantum jump probability is $\delta p(0)=0$,
i.e., no jump can occur. After integrating the system of equations by a small
step $\Delta t$ we obtain

$\displaystyle S^{z}(\Delta t)=$ $\displaystyle-1-g\Delta t+\mathcal{O}(\Delta
t^{2}),$ (75) $\displaystyle S^{0}(\Delta t)=$ $\displaystyle
1+\mathcal{O}(\Delta t^{2})$ (76)

with corresponding jump probability $\delta p(\Delta t)=-g\gamma\Delta
t/2+\mathcal{O}(\Delta t^{2})$, which is negative. For consistency, one then
has to set $\delta p(\Delta t)$ to zero.

## References

  * Binder [2005] K. Binder, Monte-carlo methods, Mathematical tools for physicists , 249 (2005).
  * Voter [2007] A. F. Voter, Radiation effects in solids, NATO Science Series (Springer, Berlin, 2007) , 1 (2007).
  * Jin  _et al._ [2016] J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio, and D. Rossini, Cluster mean-field approach to the steady-state phase diagram of dissipative spin systems, Physical Review X 6, 031011 (2016).
  * Weimer [2015] H. Weimer, Variational principle for steady states of dissipative quantum many-body systems, Physical review letters 114, 040402 (2015).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Reports on Progress in Physics 79, 096001 (2016).
  * Vidal [2004] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Physical review letters 93, 040502 (2004).
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Physical review letters 93, 207204 (2004).
  * Cui  _et al._ [2015] J. Cui, J. I. Cirac, and M. C. Bañuls, Variational matrix product operators for the steady state of dissipative quantum systems, Physical review letters 114, 220601 (2015).
  * Wootters [1987] W. K. Wootters, A wigner-function formulation of finite-state quantum mechanics, Annals of Physics 176, 1 (1987).
  * Schachenmayer  _et al._ [2015] J. Schachenmayer, A. Pikovski, and A. M. Rey, Many-body quantum spin dynamics with monte carlo trajectories on a discrete phase space, Phys. Rev. X 5, 011022 (2015).
  * Steel  _et al._ [1998] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls, and R. Graham, Dynamical quantum noise in trapped bose-einstein condensates, Phys. Rev. A 58, 4824 (1998).
  * Gardiner  _et al._ [2004] C. Gardiner, P. Zoller, and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer Science & Business Media, 2004).
  * Blakie  _et al._ [2008] P. Blakie, A. Bradley, D. M.J., B. R.J., and G. C.W., Dynamics and statistical mechanics of ultra-cold bose gases using c-field techniques, Advances in Physics 57, 363 (2008).
  * Polkovnikov [2010] A. Polkovnikov, Phase space representation of quantum dynamics, Ann. of Physics 325, 1790 (2010).
  * Zhu  _et al._ [2019] B. Zhu, A. M. Rey, and J. Schachenmayer, A generalized phase space approach for solving quantum spin dynamics, New Journal of Physics 21, 082001 (2019).
  * Perlin  _et al._ [2020] M. A. Perlin, C. Qu, and A. M. Rey, Spin squeezing with short-range spin-exchange interactions, Physical Review Letters 125, 223401 (2020).
  * Czischek  _et al._ [2018] S. Czischek, M. Gärttner, M. Oberthaler, M. Kastner, and T. Gasenzer, Quenches near criticality of the quantum ising chain—power and limitations of the discrete truncated wigner approximation, Quantum Science and Technology 4, 014006 (2018).
  * Khasseh  _et al._ [2020] R. Khasseh, A. Russomanno, M. Schmitt, M. Heyl, and R. Fazio, Discrete truncated wigner approach to dynamical phase transitions in ising models after a quantum quench, Physical Review B 102, 014303 (2020).
  * Sundar  _et al._ [2019] B. Sundar, K. C. Wang, and K. R. A. Hazzard, Analysis of continuous and discrete wigner approximations for spin dynamics, Phys. Rev. A 99, 043627 (2019).
  * Kunimi  _et al._ [2021] M. Kunimi, K. Nagao, S. Goto, and I. Danshita, Performance evaluation of the discrete truncated wigner approximation for quench dynamics of quantum spin systems with long-range interactions, Phys. Rev. Research 3, 013060 (2021).
  * Polkovnikov [2003] A. Polkovnikov, Quantum corrections to the dynamics of interacting bosons: Beyond the truncated wigner approximation, Physical Review A 68, 053604 (2003).
  * Huber  _et al._ [2022] J. Huber, A. M. Rey, and P. Rabl, Realistic simulations of spin squeezing and cooperative coupling effects in large ensembles of interacting two-level systems, Phys. Rev. A 105, 013716 (2022).
  * Huber  _et al._ [2021] J. Huber, P. Kirton, and P. Rabl, Phase-space methods for simulating the dissipative many-body dynamics of collective spin systems, SciPost Physics 10, 045 (2021).
  * Singh and Weimer [2022] V. P. Singh and H. Weimer, Driven-dissipative criticality within the discrete truncated wigner approximation, Phys. Rev. Lett. 128, 200602 (2022).
  * Brif and Mann [1999] C. Brif and A. Mann, Phase-space formulation of quantum mechanics and quantum-state reconstruction for physical systems with lie-group symmetries, Phys. Rev. A 59, 971 (1999).
  * Tilma  _et al._ [2016] T. Tilma, M. J. Everitt, J. H. Samson, W. J. Munro, and K. Nemoto, Wigner functions for arbitrary quantum systems, Phys. Rev. Lett. 117, 180401 (2016).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Reviews of modern physics 82, 2313 (2010).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, A rydberg quantum simulator, Nature Physics 6, 382 (2010).
  * Wigner [1932] E. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40, 749 (1932).
  * Fano [1957] U. Fano, Description of states in quantum mechanics by density matrix and operator techniques, Rev. Mod. Phys. 29, 74 (1957).
  * Hillery  _et al._ [1984] M. Hillery, R.F.O’Connell, M. Scully, and E. Wigner, Distribution functions in physics: Fundamentals, Phys. Rep. 106, 121 (1984).
  * Risken and Frank [1996] H. Risken and T. Frank,  _The Fokker-Planck Equation: Methods of Solution and Applications_ , Springer Series in Synergetics (Springer Berlin Heidelberg, 1996).
  * Klimov and Chumakov [2009] A. B. Klimov and S. M. Chumakov,  _A group-theoretical approach to quantum optics: Models of atom-field interactions_ (John Wiley & Sons, 2009).
  * Gardiner and Zoller [2010] C. Gardiner and P. Zoller,  _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ , Springer Series in Synergetics (Springer Berlin Heidelberg, 2010).
  * Žunkovič [2015] B. Žunkovič, Continuous phase-space methods on discrete phase spaces, EPL (Europhysics Letters) 112, 10003 (2015).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually-controlled rydberg atoms, Nature Phys. 16, 132 (2020).
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia, Journal of Open Research Software 5 (2017).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review 59, 65 (2017).
  * Pucci  _et al._ [2016] L. Pucci, A. Roy, and M. Kastner, Simulation of quantum spin dynamics by phase space sampling of bogoliubov-born-green-kirkwood-yvon trajectories, Phys. Rev. B 93, 174302 (2016).
  * Höning  _et al._ [2013] M. Höning, D. Muth, D. Petrosyan, and M. Fleischhauer, Steady-state crystallization of rydberg excitations in an optically driven lattice gas, Phys. Rev. A 87, 023401 (2013).