Variational Dynamics of Open Quantum Spin Systems in Phase Space

Abstract

We introduce a variational method for simulating the dynamics of interacting open quantum spin systems. The method is based on the spin phase-space representation and variationally targets the Husimi-$Q$ function with an ansatz based on a multi-dimensional mixture of spin-coherent states. Crucially, the mixture coefficients are allowed to take negative values, enabling the faithful capture of quantum correlations beyond semiclassical descriptions. The resulting equations of motion are derived from the Dirac-Frenkel variational principle and can be evaluated efficiently without resorting to Monte Carlo sampling by exploiting the analytical structure of the ansatz. As a first application, we demonstrate that this approach accurately captures both the full quantum dynamics and the non-equilibrium steady states of the transverse-field quantum Ising model, in excellent agreement with exact diagonalization. Furthermore, we show that the method scales efficiently to large two-dimensional lattices, a regime that remains challenging for other techniques.

Figures (4)

• Figure 1: Schematic of the variational phase-space approach introduced in this work. (a) A driven-dissipative spin lattice with $N_s$ sites governed by the Lindblad master equation $\partial_t \hat{\rho} = \mathcal{L}\hat{\rho}$, with nearest-neighbour $ZZ$ coupling of strength $V$, transverse-field $g$, and local decay with rate $\gamma$. (b) The many-body density matrix is encoded in phase space via the Husimi-$Q$ function, which is parameterized as a variational mixture of $N_c$ spin-coherent-state products: $Q(\bm{\Omega};\bm{\theta}) = \sum_{k=1}^{N_c} c_k \prod_{i=1}^{N_s} q(\bm{n}_i;\bm{m}_{ki})$, where each single-site function $q({\boldsymbol{n}}_i;\bm{m}_{ki})$ [cf. Eq. \ref{['eq:vcms']}] lives on the Bloch sphere for site $i$. (c) The phase-space Liouvillian $\mathcal{L}_Q$ generates the dynamics of $Q(\bm{\Omega};\bm{\theta})$; projecting onto the variational manifold via the Dirac-Frenkel principle yields the equations of motion $T\dot{\bm{\theta}} = \bm{F}$, where the quantum geometric tensor $T$ and the force vector $\bm{F}$ are evaluated through analytical integrals [cf. Eq. \ref{['eq:S-F']}] and automatic differentiation, without the need of Monte Carlo sampling.
• Figure 2: Variational steady-state results for the observables $\langle \hat{\sigma}_x\rangle$ and $\langle \hat{\sigma}_y\rangle$ for a $1$D $16\times 1$ dissipative transverse-field Ising model. The parameter $g$ represents the external magnetic field along $X$ while $\gamma$ is the loss rate associated to the jump operator $\hat{\Gamma}_i = \sqrt{\gamma} \hat{\sigma}_i^-$. The v-MCS is compared to neural-network steady-state results with Transformer and Convolutional Neural Networks (CNN) architectures (data taken respectively from Luo2022 and Mellak2024, see legend). Exact results are obtained via Monte Carlo wave function simulations. The v-MCS ansatz at convergence has $N_c =16$ total components per spin for a total of $784$ parameters. The initial condition for the dynamics is $\prod_{i=1}^{N_s}\ket{\uparrow}_x$. The $ZZ$ coupling strength is $V /\gamma = 2$.
• Figure 3: Benchmark results in 2D spin lattices. (a) Variational real-time dynamics (solid lines) of a $3\times 3$ dissipative transverse-field Ising model for $g/\gamma = 1, 2, 3$, compared with the exact solution (dashed lines). The v-MCS ansatz uses $N_c = 10$ components per spin ($280$ variational parameters). The initial state is a small perturbation of $\prod_{i=1}^{N_s}\ket{\uparrow}_x$. (b) Steady-state expectation values $\langle \hat{\sigma}^\alpha \rangle_\mathrm{ss}$ for a $4\times 4$ lattice as a function of $g/\gamma$. The v-MCS results ($N_c = 16$, $784$ parameters) are compared with exact results (see legend). In both panels, $V/\gamma = 2$.
• Figure 4: Variational dynamics of $\langle\hat{\sigma}^\alpha\rangle$, with $\alpha \in \{x,y,z\}$ for the $8\times 8$ dissipative transverse-field Ising model, computed with $N_c \in \{2, 4, 6\}$ coherent-state components (respectively with $386,772$ and $1158$ variational parameters). System parameters: $V/\gamma = 2$, $g/\gamma = 2$. The initial state is $\prod_{i=1}^{N_s}\ket{\uparrow}_x$. Both transient dynamics and steady-state values exhibit systematically controlled convergence as $N_c$ increases.