Many-Body Quantum Spin Dynamics with Monte Carlo Trajectories on a Discrete Phase Space

Interacting spin systems are of fundamental relevance in different areas of physics, as well as in quantum information science, and biology. These spin models represent the simplest, yet not fully understood, manifestation of quantum many-body systems. An important outstanding problem is the efficient numerical computation of dynamics in large spin systems. Here we propose a new semiclassical method to study many-body spin dynamics in generic spin lattice models. The method is based on a discrete Monte Carlo sampling in phase-space in the framework of the so-called truncated Wigner approximation. Comparisons with analytical and numerically exact calculations demonstrate the power of the technique. They show that it correctly reproduces the dynamics of one- and two-point correlations and spin squeezing at short times, thus capturing entanglement. Our results open the possibility to study the quantum dynamics accessible to recent experiments in regimes where other numerical methods are inapplicable.


I. INTRODUCTION
The study of spin models with long-range interactions is of fundamental relevance in different areas of physics ranging from condensed matter to atomic and molecular physics, quantum optics, quantum information science and even biology. Recent experimental progress is allowing us to access key aspects such as the build-up of long-range correlations and entanglement, and the propagation of information during non-equilibrium dynamics [1][2][3]. However, many questions remain open, partly due to the absence of appropriate theoretical tools to calculate the time-evolution in complex quantum systems. For example, numerical time-dependent density matrix renormalization group (t-DMRG) methods [4][5][6] become inefficient in higher dimensional systems. Traditional perturbative and Keldysh techniques [7], as well as kinetic theories and phase space methods [8][9][10] are limited to weakly interacting, close to equilibrium or short time situations. Cluster expansions [11] are limited to dilute samples with moderately short range interactions. Consequently, of immediate relevance is the development of new numerical techniques capable of overcoming many of these shortcomings. In this work, we advance towards this direction by introducing an improved phase-space method that we refer to as the discrete Truncated Wigner Approximation (dTWA).
The dTWA is based on a sampling of a discrete Wigner function. We show that it greatly outperforms other currently used phase-space methods, such as the Truncated Wigner Approximation (TWA). Standard phase space methods replace the quantum-mechanical time-evolution by a semi-classical evolution via classical trajectories. The quantum uncertainty in the initial state is accounted for by an average over different initial conditions [8,9] determined by a continuous Wigner function. In contrast, the use of discrete Wigner functions enables us to quantitatively access dynamics in generic spin lattice models, including oscillations and revivals of single particle observables and correlation functions not captured by the TWA. As a particular example, we show that dTWA can capture the build-up of spin squeezing [12,13] in spinmodels with finite-range interactions.
This paper is organized as follows: In Sec. II the traditional TWA is reviewed, and in Sec. III our dTWA technique is introduced. Sec. IV contains a benchmark of the drastic improvement provided by the dTWA by comparing dynamics of single-spin observables, correlation functions and spin-squeezing for a model with Ising and XY interactions. Sec. V provides a conclusion and an outlook.

II. TRUNCATED WIGNER APPROXIMATION (TWA)
The mapping between the Hilbert space of a quantum system and its corresponding classical phase space can be accomplished through the so-called "phase-point operators"Â. Those relate the density matrix of the quantum systemρ to a classical quasi-probability distribution (generally non-positive) known as Wigner function W [28]. The phase-point operators can be written in terms of phase space variables p and q as (we set = 1 in this paper) [28,29]: with the corresponding Wigner function. In general however it is not possible to accomplish this exactly. The TWA [8,9] approximates the dynamics by taking only first order quantum fluctuations into account. This is done by assuming both that the Wigner function is conserved on classical trajectories fixing it to its initial state, W (p, q) → W (p 0 , q 0 ), and that the Weyl symbols follow a classical evolution, obtained by solving Hamilton's equations of motion, O W (p(t), q(t)) → O W (p cl (t), q cl (t)): Generalizations of this continuous formulation to N spin-1/2 particles have been developed, e.g. by means of a spin-coherent state representation, valid up to 1/N corrections. Generically, the Wigner function is approximated by a smooth positive Gaussian-like distribution in the collective spin variables. For example if all spins are pointing along the z-axis, the Wigner function can be approximately written as [10] W (S ⊥ , S z ) ≈ 1/(πS)e S 2 ⊥ /S δ(S z −S), with S = N/2, S ⊥ = (S 2 x +S 2 y ) 1/2 , and S z the transverse and longitudinal spin component of the collective classical spin, respectively. This Wigner function has a clear interpretation: If each spin initially points along the z direction, because of the Heisenberg uncertainty principle the transverse spin components must fluctuate as S 2 ⊥ ∼ S. This Gaussian TWA generally captures aspects of the quantum spin dynamics at short times but lacks important ones such as the revivals, ubiquitous in discrete quantum systems, and it is mainly limited to deal with the collective dynamics [9]. For interactions with spatial structure the dynamics quickly takes the system out of the collective-spin Hilbert-space and the TWA breaks down. In the continuous phase-space picture ways to overcome these shortcomings have been proposed using hidden variables [30] or more complex representations of the Wigner function [9,31,32]. Here we propose and test a completely different ansatz by employing discrete phase-spaces [29].

III. DISCRETE TRUNCATED WIGNER APPROXIMATION (DTWA)
For systems with discrete degrees of freedom, it is possible to define a "discrete phase-space" in many ways (see [29] and references therein). We will use four phasepoints to describe a single spin-1/2, which we denote as α ≡ (q, p) ∈ {(0, 0), (0, 1), (1, 0), (1, 1)} as introduced by Wootters [29,33]. For each phase space point one can define a phase-point operatorÂ α , the Weyl symbols O W α = tr(ÔÂ α )/2 (or inverselyÔ = αÂ α O W α ), and the Wigner function w α which is the Weyl symbol for the density matrix. A set of phase-point operators is The quantum physics of a spin-1/2 particle can be fully described by a discrete four-point Wigner quasi-probability distribution, w (p,q) . The probability for a spin to point along the ±z, ±y, and ±x direction (p x,y,z ±1 ) is given by the sum over the vertical, diagonal, and horizontal lines in phase space, respectively [29]. (b) Example of a spin pointing along +z. (c) The concept of the dTWA for a system of N spins is to i) randomly sample the phase points for each spin according to w [n] p,q ; ii) evolve the phase points according to classical equations; and iii) evaluate observables in phase space from the statistical mixture of nt trajectories.
given by [29]: Hereσ = (σ x ,σ y ,σ z ) are the Pauli matrices. In a many-body system with N spin-1/2 particles the discrete phase space has 4 N points we denote as α = {α 1 , α 2 , . . . α N }. Analogously to Eq. (1) on this discrete phase-space we can now formulate the dTWA as: where w α (0) is the initial Wigner function on the discrete many-body phase space and O W,cl. α (t) the classically evolved Weyl symbol (see Appendix A for more on the classical equations of motion). Numerically, we can solve Eq. (3) by statistically choosing [according to w α (0)] a large number n t of random initial spin-configurations. Each "trajectory" is evolved in time according to the classical equations of motion and the expectation value in Eq. (3) is estimated via statistical averaging (error ∼ 1/ √ n t ). We find that the required n t does not depend on the system size, but rather on the observable under consideration (see Appendix C).
As an example of how to construct the initial Wigner function, let us again consider an initial state with all spins pointing along the z direction. The initial density matrix factorizes,ρ(0) = N i=1℘ [i] (ẑ) (the superscript [i] denotes the Hilbert/phase space for spin i), and thus αi . Here, w (0,1) = 1/2, and w (1,1) = 0 for every spin i (cf. Fig. 1 for an illustration). Note that for this initial state all quasi-probabilities are positive, which is important for the numerical sampling. The sum along each "phase-space line" of the discrete Wigner function can be associated with the probability of a measurement outcome [29], similarly to the continuous variable case. As shown in Fig. 1, here w (0,1) = 1/2 means that the probability to find an individual spin pointing along the +z-direction is one (sum over the first row), while the probabilities to find it along +x or −x (sum over each column) are 50% each; as well as for +y or −y (sum over the diagonals). Note that the discrete sampling properly accounts for the quantum mechanical probability distribution of the x, y, z spin components of a qubit in the sense that all moments are reproduced correctly: (σ x,y ) k = (1 + (−1) k )/2, (σ z ) k = 1, with k a positive integer.

From Eq. (3) and the classical equations of motion (Appendix
αi (t) , which resembles the time-evolution of an expectation value computed from a separable density matrix. However, the dTWA is able to capture the build-up of spin-squeezing (an entanglement witness [34]). This is not a contradiction, because the phase-point operatorsÂ α are in fact not density matrices. While Tr(Â α ) = 1, it is possible to satisfy Tr(Â 2 α ) > 1; these conditions might be interpreted as a statistical mixture with negative probabilities.

IV. RAMSEY DYNAMICS
To demonstrate the accuracy of the dTWA we consider a system of N two-level systems arranged on a lattice with M sites with dynamics governed by the Hamiltonian:Ĥ We consider an initial product state in which all spins are aligned along the x axis. The interactions are assumed to decay as a function of the distance with a decay exponent α, and allowed to be spatially inhomogeneous: e.g. J Here, r ij is the vector connecting spins on sites i and j and θ is the angle it makes with the quantization axis (chosen along z). We discuss two specific cases, Ising (J ⊥ ij = 0) and XY (J z ij = 0) interactions. The Ising limit is naturally realized in experiments with ion traps (both in 1D [19,20] or 2D [21] geometries) and Rydberg atoms in 2D [17]; the XY limit dynamics has been realized in polar molecules in optical lattices [11,14], and magnetic atoms [35], Förster resonances in Rydberg atoms [36], and as an effective Hamiltonian in trapped ions with a superimposed large transverse field [20].

A. Ising interactions
In this limit, exact analytical expressions for the dynamics exist [37][38][39] and can be used to benchmark the dTWA. Dynamics of observables involving the collective spin S = n σ n for a system with M = N = 100 sites in a 1D chain (oriented along x) are shown in Fig. 2. We calculate the time-evolution of S x as well as the collective correlation functions ∆S x = S 2 x − S x 2 and Re S y S z . We consider the case of all-to-all (α = 0) and shortranged dipolar interactions (α = 3). In the all-to-all case S x shows revivals at times that are multiples of π/2J. In contrast to the traditional Gaussian TWA which only captures the initial decay of the contrast and misses the revivals, the dTWA fully reproduces the exact dynamics. The α = 3 case exhibits an oscillatory dynamics, perfectly accounted for by the dTWA solution, but not captured by the traditional TWA. The dTWA calculations for Re S y S z also show perfect agreement. For the correlation ∆S x deviations are visible, however, the oscillatory dynamics in ∆S x is still better reproduced in dTWA than in the traditional TWA.
We understand the agreement between the Ising solution and the dTWA analytically. For a particular spin n, the exact solution for the time-evolution is σ x n (t) =  N = 100). The chain is aligned along the x direction (θ = π/2). The dTWA gives exact results for short times and good estimates for the achievable ξ. (b) Achievable ξ as a function of the filling fraction (averaged over 1000 random configurations). Exact diagonalization (ED) for N < 20, tDMRG otherwise. Inset: rapid increase of the entanglement entropy for non-integer fillingn < 1 (single configuration), rendering t-DMRG calculations inefficient in this regime.
the z-component of the spin to be randomly distributed as s z i (t = 0) = ±1, we can identify s z i (t = 0) with the parameters m i in the exact solution and find that the dTWA becomes equivalent to the exact solution for large n t and arbitrary numbers of spins N . Note that in addition the traditional TWA approach is only valid in the large-N limit. The comparisons with the exact solution, extended to correlation functions, confirm the excellent agreement of Re S y S z and show the origin of the discrepancies in other two-point correlations (see Appendix B and Fig. 2).

B. XY model
Here we study systems with varying filling fractions n = N/M , a situation relevant for recent polar molecule experiments [14]. We focus our attention on the timeevolution of spin-squeezing, which is a signature of quantum correlations and entanglement [34] and is a resource for enhanced sensitivity in quantum metrology [40]. The spin squeezing parameter is ξ ≡ √ N min n ⊥ (∆S ⊥ )/|S| where S is the total collective spin, S ⊥ = S · n ⊥ , and the minimum is taken over all unit vectors n ⊥ (directions) perpendicular to the vector S. On the Bloch sphere a squeezed state with ξ < 1 shows an elliptical profile of the spin noise distribution [13].
For the XY model no exact solution exists for generic spin systems, however in 1D we can use t-DMRG to calculate the exact dynamics at short times. This technique works as long as the bi-partite entanglement in the system, quantified by the half-chain von Neumann entropy S vN (ρ L ) = −tr(ρ L log 2 ρ L ) (where ρ L is the reduced density matrix of half the spin-chain), remains small [41][42][43][44]. Surprisingly, we find that for intermediate filling fractions, due to the inhomogeneity in J ⊥ ij , the entropy S vN grows much more rapidly than forn = 1 (see inset in Fig. 3b). Thus with reasonable computational resources exact results could only be calculated for (N ≤ 32) and (N ≤ 100) in a system with M = 100 sites. Our results for a 1D chain of spins along the x direction are shown in Fig. 3. In Fig. 3a we compare the evolution of ξ from an exact t-DMRG calculation with the dTWA result forn = 1. We find that dTWA agrees on short times and captures almost all of the spin squeezing that is created in the time-evolution. In Fig. 3b we plot the maximally achievable spin squeezing as function of the filling fraction. The dTWA interpolates over the whole range of filling fractions connecting the exactly tractable limits. It is interesting to note that while squeezing implies entanglement (non-separability) [34], the type of squeezing we consider here is apparently independent of bipartite entanglement, in the sense that for smaller S vN (ρ L ) we can have large squeezing and vice versa.

C. Varying range of interactions and dimension
We now check the dependence of the validity of the dTWA on the range of interactions and the dimensionality of the system, focusing on the XY model. The results are summarized in Fig. 4. In Fig. 4a we analyze the time-evolution of ξ for a M = N = 20 spin system in 1D for various α and for fixed α and different dimensions. As expected we find that longer-ranged interactions lead lo larger amounts of spin-squeezing. Remarkably, for longer-ranged interactions the increased squeezing is better captured within the dTWA (also when analyzing the relative error). This can be understood in the limit of all-to-all interactions (α = 0), where due to the conservation of the total spin the XY model becomes equivalent to the Ising model. In that case we found (Fig. 2) that the dTWA is almost exact. For higher dimensional systems, we find that the agreement becomes generally better. However due to the inhomogeneity of the coupling constants, the 3D case is difficult to access since the spinsqueezing becomes very small. In Fig. 4b we find that for a small 2D system with α = 1, the agreement is excellent.

V. CONCLUSION & OUTLOOK
The discrete Truncated Wigner Approximation (dTWA) opens the possibility of computing the dynamics of large spin systems in regimes where at the present time there is no other theoretical tool at hand (see Fig. 4b). This is possible since the computational time for a solution of the mean-field equations only scales polynomially in time. Furthermore, the Monte-Carlo sampling can be parallelized, and the number of required samples for statistical convergence depends only on the observable (Appendix C). In the future it will be interesting to try to combine higher order corrections to the TWA [9,30] with the idea of discrete phase spaces. This could lead to methods that are able to also capture dynamics for longer time-scales.

APPENDIX B: CORRELATION FUNCTIONS IN THE ISING MODEL
The time-evolution for the Ising Hamiltonian (6) can be solved exactly (see references in the main text). The time-evolution of the local operator at site k [σ ± k = (σ x k ± iσ y k )/2] can be calculated as: Here each m i takes the values −1 and +1, and the sum runs over all 2 N possible combinations. Comparing with Eq. (10) one sees that dTWA gives the exact time evolution in this case, when the sum is approximated via a random sampling of s z taking the values +1, −1 (see discussion in the main text). The same calculation is possible for correlations. For example, between particle i and j (i < j): with 1 , 2 ∈ {+1, −1}. Note that the two sums with indices m i and m j are missing. This is to be compared with the dTWA which estimates the same quantity (in the limit n t → ∞) as One sees that the only difference between (12) and (13) are the two additional sums with s z i , s z j ∈ {−1, +1} in (13). This gives rise to an error of Since for example S 2 x /N 2 = 1 N + 1 c. , this correlation contains an error. Analogously, it is straightforward to see that due to the conservation of the z-component of the spin there is no error made when calculating Re S z S y = i =j 2Re σ z j σ − i with dTWA.