A generalized phase space approach for solving quantum spin dynamics

Numerical techniques to efficiently model out-of-equilibrium dynamics in interacting quantum many-body systems are key for advancing our capability to harness and understand complex quantum matter. Here we propose a new numerical approach which we refer to as GDTWA. It is based on a discrete semi-classical phase-space sampling and allows to investigate quantum dynamics in lattice spin systems with arbitrary $S\geq 1/2$. We show that the GDTWA can accurately simulate dynamics of large ensembles in arbitrary dimensions. We apply it for $S>1/2$ spin-models with dipolar long-range interactions, a scenario arising in recent experiments with magnetic atoms. We show that the method can capture beyond mean-field effects, not only at short times, but it also correctly reproduces long time quantum-thermalization dynamics. We benchmark the method with exact diagonalization in small systems, with perturbation theory for short times, and with analytical predictions made for closed system which feature quantum-thermalization at long times. By computing the Renyi entropy, currently an experimentally accessible quantifier of entanglement, we reveal that large $S$ systems can feature larger entanglement than corresponding $S=1/2$ systems. Our analyses demonstrate that the GDTWA can be a powerful tool for modeling complex spin dynamics in regimes where other state-of-the art numerical methods fail.


I. INTRODUCTION
In the past years, rapid developments of various experimental platforms have made it possible to observe outof-equilibrium dynamics of large isolated quantum manybody models in controlled environments [1][2][3][4]. Naturally, this also leads to a high demand for numerical methods capable of simulating such dynamics. Computations for large system sizes beyond a classical mean-field picture are a challenging task due to the complexity of the full quantum problem. Consequently, most recently developed methods for large systems are limited to either low dimensional geometries (e.g. by making use of matrix product states/tensor networks [5][6][7][8][9][10]), particular ansatz wave-functions (e.g. in combination with variational Monte-Carlo evolution [11][12][13][14]), or dilute systems (e.g. by making use of a clusterization [15,16]).
In particular, models of coupled spin-particles with long-range interactions have become a topic of intensive research because of important experimental progress. While models with spin S = 1/2 have been implemented with many different setups, e.g. using polar molecules [15,17], Rydberg atoms [16,[18][19][20][21][22][23], trapped ions [24][25][26] and cavity QED systems [27,28], recently also models with larger spins S > 1/2 have become a research focus in particular for experiments with magnetic atoms [29][30][31][32][33][34]. The large spin degrees of freedom in S > 1/2 systems poses a much more stringent requirement for numerical treatment compared to S = 1/2 systems. The full Hilbert space involves (2S + 1) N states for N particles, which renders exact diagonalization (ED) impossible already for small N . In addition, many current experiments are performed in 2D or 3D, thus also disabling otherwise very successful matrix product state techniques for long-range interacting systems in 1D [35][36][37]. This calls for new theoretical tools capable of accounting for the many-body nature of these models as well as intrinsic quantum correlations.
In this paper, we present an efficient numerical approach based on a semi-classical phase space method that can be applied to large 2D/3D lattice systems with arbitrary spin S. The method is based on the well known truncated Wigner approximation (TWA) [38][39][40][41] adapted to spin-models. The general TWA idea relies on a sampling of the quantum fluctuations of the initial state from a Wigner function, and an evolution of the samples along classical trajectories. Importantly, in contrast to previous approaches, here we introduce an enlarged phase space representing not only 3 spin-variables, but all (2S + 1) 2 density matrix elements for each spin [29,34,42,43]. Furthermore, we can use a sampling from discrete quasi-probability distributions, which are exact and positive for most experimentally relevant initial states. Our method thus allows us to study the TWA time-evolution not only of spin operators, but the full spin-density matrix and thus allows us to extract experimentally relevant time-dependent observables such as spin-state populations, and fundamentally relevant quantities such as entanglement. In the limit of S = 1/2 our generelized discrete TWA approach (GDTWA), reduces to the previously proposed discrete TWA method (DTWA) [44], which has been remarkably succesful in predicting S = 1/2 model dynamics [16,[45][46][47][48][49]. This paper is organized as follows: In Sec. II, after a brief review of the TWA for continuous Wigner functions, we introduce the GDTWA. In Sec. III, we test its validity by a comparison with ED. We focus on the evolution of Zeeman level populations induced by dipolar interactions (relevant to experiments using both Cr [29,[31][32][33] and Er atoms [34]) and the evolution of entanglement. In Sec. IV, we apply the GDTWA to investigate various aspects of spin dynamics: the spreading of population in a synthetic dimension and the underlying approach to thermalization, as well as the build-up of quantum entanglement. Finally, in Sec. V we conclude and discuss applications for other systems as an outlook.

II. METHOD
A. Phase space sampling for quantum spin systems The Hamiltonian for a system consisting of N spin-S particles coupling to each other via two-body interactions can be written aŝ with α, β = x, y, z, and i, j = 1, 2, ...N . The first term governs local fields and the second term inter-spin interactions. The dynamics of the three components of the spin-operator of the i particle, Ŝ i α=x,y,z , can be obtained via the Heisenberg equations of motion (we set ≡ 1 throughout this paper): These equations generally depend on high order interspin correlations such as Ŝ i αŜ j β , Ŝ i αŜ j βŜ k γ , . . . , for i = j = k and in practice it is impossible to solve them exactly for large N . In a mean-field ansatz one assumes that such spin correlations factorize and can be written in terms of single particle observables, Ŝ i α , e.g., Ŝ i αŜ j β ≈ Ŝ i α Ŝ j β . Then Eqs. (2) turn into N coupled closed equations which can be easily solved numerically, and which correspond to equations for classical spin-variables 1 . Importantly, the factorization neglects entanglement between the spins, as the total state of the system is forced to remain a product-state. The missing quantum correlations are in many cases crucial. One general approach to account for some of those correlations is to retain higher order correlations using e.g. BBGKY 1 Note that more general models than Eq. (1) for S > 1/2 can also depend on intra-spin terms such as quadratic fields ∝ (Ŝ i z ) 2 . In this case also the mean-field equations are not necessarily closed unless also intra-spin correlations are assumed to factorize. Those cases are not accounted for in the approach described here, but we will properly include them in our GDTWA method below.
hierarchies [46]. However, then already retaining secondorder correlations increases the complexity of Eqs. (2) to ∼ N 2 and furthermore higher order corrections can typically lead to numerical instabilities. An alternative approach to simulate quantum correlations while retaining a complexity ∼ N is to resort to a phase space description of the quantum system [40,41] as reviewed in the following.
In the phase space approach, quantum operatorsÔ are mapped to functions in the phase space of classical variables, so-called Weyl symbols. Taking as an example a particle moving in 1D, the phase space variables are given by position and momentum x and p, respectively, and the Weyl symbol is denoted as O W (x, p). The expectation value of any operatorÔ can then be computed as an integral in phase space: Here, W (x, p) is the Wigner function corresponding to the Weyl symbol of the density matrixρ: W (x, p) = 1/(2π) ds x + s/2|ρ |x − s/2 exp(−ips), and represents a quasiprobability distribution for points in phase space. In the truncated Wigner approximation (TWA), the time evolution of Ô is approximated from the initial Wigner function and the classical evolution of the phase space variables: Here, x cl (t) and p cl (t) are classical trajectories of x and p, respectively, with x cl (0) = x 0 and p cl (0) = p 0 . Then, compared to a fully classical evolution, in the TWA dynamics quantum fluctuations are taken into account to the lowest order, in the sense of keeping only terms linear in in the equations of motions for the Weyl symbols [40,41]. Observables (hermitian operators) give rise to real Weyl symbols and thus correspond to symmetrized sums of the position and momentum operator. For example, the Weyl symbol O W (x cl (t), p cl (t)) = x 2 cl (t)p cl (t) corresponds to the observableÔ = (x 2p + 2xpx +px 2 )/2 [41].
For a system of N coupled spin-S particles, a natural way to formulate the TWA dynamics consists of replacing {x,p} → {Ŝ i α }, and using the 3N spin-variables {S i α } (with α = x, y, z and i = 1, 2, . . . N ), corresponding to the Weyl symbols of {Ŝ i α }, as phase space variables. The trajectories S i α,cl (t) correspond to those obtained from the mean-field approximation of Eq. (2) with S i α,cl ≡ Ŝ i α . For particular initial states, the Wigner function can be easily found. Taking for example a product state iρ i , where each spin i is described by a spin coherent state, then for large S the Wigner function can be well approximated by a Gaussian. In the case of a state initially polarized along the z direction, the Wigner function factorizes, and for each spin takes the simple where δ(·) denotes the delta function. For this state the phase space value of S i z is determined as S, while the values of S i x (S i y ) fluctuate according to a Gaussian distribution with zero mean and a variance ∆(S i x ) 2 = ∆(S i y ) 2 = S/2 (see Fig. 1). This width of the distribution of the classical variables reflects the uncertainty relation intrinsic to the quantum mechanical operatorsŜ i x andŜ i y . For large S 1, ∆S i x /S → 0, suggesting vanishing effects from quantum fluctuations.
For such initial states, the TWA evolution now reads where the Weyl symbol corresponds again to symmetrized observables. For example, for an inter-spin correlation such asÔ =Ŝ n xŜ m x at distinct sites n = m, It is important to note that, in contrast to a mean-field simulation, a TWA evolution can lead to inter-spin quantum correlations. For example, due to the time-evolution according to (generally) non-linear classical equations, a spincomponent for a single spin m can now depend on the initial conditions of all other spin-variables via some function, e.g. S m x,cl (t) ≡ F m x ({S j β,0 }). Therefore, for observables computed in the TWA, such asŜ m xŜ n x with m = n, There are some important shortcomings of the Gaussian TWA approach with 3 spin-variables introduced in the previous section. Most importantly, the scheme is tailored to problems relying on the 3 spin operators, both in the Hamiltonian and for measurements. Most experimental scenarios go beyond this limitation, by including e.g. quadratic fields ∝ (Ŝ i z ) 2 and measurements of Zeeman state populations [29,34]. Furthermore, the approach is insufficient for a full state-tomography in large S systems (see B). More generally, we would like to adapt the method to arbitrary many-body models of coupled discrete local Hilbert spaces.
To go beyond the limitations we proceed by enlarging our phase space from 3 spin-variables to elements of higher-dimensional generalizations of Bloch vectors [50]. This can be accomplished by noticing that for a spin-S atom with N = 2S + 1 spin states, its density matrixρ i consists of D = N × N elements. Correspondingly, we can define D hermitian operators, Λ µ (hats are dropped for clarity of notation), using the generalized Gell-Mann matrices (GGM) and the identity matrix 1 2 : For each spin i, these matrices are orthonormal, Tr[Λ i µ Λ i ν ] = δ µ,ν , and constitute a complete local basis. Hence any single-spin operator can be represented viâ and µ = 1, 2, ..., D. Consequently, more general Hamiltonians than Eq. (1), with one-and two-body terms can now be represented aŝ with µ, ν = 1, 2, . . . , D and i, j = 1, 2, . . . , N . The evolution of expectation values of GGMs can now be computed as For each GGM Λ i µ we define a corresponding real phase space variable, λ i µ . As in the ordinary TWA we will assume that the dynamics can be determined from statistical averages over trajectories of the phase space variables only. Thus in analogy to the case of 3 spin-variables discussed after Eq. (2), our classical equations of motions follow from assuming a factorization between different sites Λ i µ Λ j ν . . . Λ k σ = Λ i µ Λ j ν . . . Λ k σ for any nonequal site-index i, j, . . . , k in Eq. (14), and by replacing λ i µ (t) ≡ Λ i µ . To define a phase space probability distribution for the inital state, we decompose each Λ i µ via its eigenvectors |a i µ with corresponding eigenvalues a i µ , Λ i µ = a i µ a i µ |a i µ a i µ |. As in spin-1/2 systems, where the Pauli matrices σ x,y,z can be measured to be ±1 in a projective measurement, the a i µ correspond to possible measurement outcomes of Λ i µ . We will focus on the experimentally relevant case of initial product statesρ = iρ i . Then, for eachρ i we can define a corresponding probability distribution for our phase space variables, which we limit to a discrete set of values The overall distribution factorizes for different variables on the same site and between sites, such that the probability for a configuration of all λ i µ being a certain combination of the eigenvalues a i µ for i = 1, 2, . . . , N and µ = 1, 2, ..., D − 1 is given by From the distribution (15), we can compute the expectation value of any observable on a fixed site m via the expansion (12), Here we defined the notation · ≡ , denoting the statistical average over any combination of eigenvalues for the λ i µ inside the brackets. Note that for a Gaussian distribution for three spin-variables it is not possible to exactly represent any observable on the Hilbert space of a single-spin, since the three spinoperators on a site m,Ŝ m x,y,z and 1 do not form a complete local basis for any operator. Thus, for example the observable (Ŝ m z ) 4 cannot be expanded as linear superposition of the matrices {Ŝ m x,y,z , 1} and has to be evaluated as 4-th order moment from the probability distribution. As a consequence, for a spin coherent state polarized along x, the Gaussian distribution does not reproduce (Ŝ m z ) 4 correctly (see B).
For our GDTWA scheme, we use a discrete configuration selected from the distribution (15) as initial condition for a classical evolution of the phase-space variables λ i µ (0) ∈ {a i µ } for all µ and i. This configuration is then numerically evolved according to the equations of motion that we derive from the factorization of Eq. (14). Those equations are in general non-linear and can be written in the form dλ(t)/dt = M(λ(t))λ(t). Here λ is a vector of all N × (D − 1) phase-space variables λ i µ , and the matrix M(λ(t)) depends in general on the configuration at time t. We solve those equations numerically for the selected initial condition and obtain a mean-field trajectory, λ m ν (t). We repeat this procedure and average over the possible initial configurations. We thus obtain e.g. an approximation for the expectation value of a GGM Λ m ν at time t as Λ m If the problem is linear, i.e. if M(λ(t)) = M is timeindependent, e.g. if there are only single-spin terms in Eq. (14), then the time-dependent solution in our scheme is given by Λ m (16), the time-evolution scheme is exact. If the problem is non-linear, which is generally the case in presence of inter-spin interaction terms in Eq. (14), then the statistical average λ m ν (t) will also depend on intra-spin correlations in the initial state, also within a single spin, e.g. λ i µ (0)λ i ν (0) , and higher orders. In general, those correlations do not coincide with the exact correlations obtained from the true quantum state. For example, it is straightforward to see that per definition of our discrete probability distribution, all correlations factorize between different GGMs on the same site, e.g.
This is not necessarily true for the quantum mechanical intra-spin correlations of a generic single-spin state. For example, for some stateŝ where we have to use a symmetrization since the GGMs do not generally commute.
In A we show, that for "diagonal" product states all initial intra-spin correlations can be perfectly reproduced by our distribution (15). We define those states aŝ ρ = i |α 0 i α 0 | where α 0 = 1, 2, . . . , N , and |α 0 i is any of the single-spin eigenstates of the matrices defined in Eqs. (10). It is important to mention that any product of non-diagonal pure states can be brought into the formρ = i |α 0 i α 0 |, via local unitary transformations. Therefore, the diagonal assumption is not an actual restriction and in our simulation scheme we can exactly describe generic products of initial pure states. This can be achieved by applying the same unitary transformations to the equations of motion (see A). Note that, alternatively, one may also use a Gaussian approximation for the distribution of the initial λ i µ [42,43], i.e. for all D − 1 generalized Bloch sphere variables. However, we point out that in contrast to the discrete distribution given by Eq. (15), the Gaussian sampling not only does not reproduce all initial intra-spin correlations correctly, but also can lead to worse longer-time predictions. For example, the dynamics of collapses and revivals has been only observed with the exact discrete sampling [51].
To compute the time-dependent expectation value of an arbitrary multi-spin observable,Ô, we expand and approximate This equation is a discrete analog of the TWA method from Eq. (6) for a phase space extended to a highdimensional generalized Bloch sphere. Approximating Eq. (18) with a finite number of samples, n s , from the discrete probability distribution is what we denote as generalized discrete truncated Wigner approximation (GDTWA). The key idea underlying the method is also illustrated in Fig. 1.
As a clarifying example, let's also explicitly provide a formula for computing the evolution of a single-spin observable,Ô m . After selecting a single random configuration for all phase space variables {λ i µ } according to , we denote the subsequent classical evolution of the variables at site m for this sample "s" as λ m, [s] µ,cl (t). Then, the evolution ofÔ m is computed as Note that in contrast to a continuous Wigner function approach, here in principle only a finite number of numerical samples n s is needed to compute Eq. (18). In practice, however, this number increases exponentially with the N . The n s needed to obtain converged results depends on the system size and the observable. In practice for typical realistic problems (with hundreds of spins as used in Sec. IV) we find a number on the order of a few 10 4 samples sufficient. Note that we find that for collective observables the number of samples necessary for convergence typically decreases with N .
To further explain the main idea behind this approach and its capabilities, we first consider as an example an array of spin S = 1/2 particles. In this case, the three nontrivial Λ µ for each spin are proportional to standard Pauli matrices, Λ x,y,z = σ x,y,z / √ 2, each of which has eigenvalues ±1/ √ 2, and eigenstates described by states aligned (anti-aligned) along the corresponding direction: |↑ x,y,z (↓ x,y,z ) . For a spin initially polarized as |↑ z , the initial values of λ x in phase space take its two eigenvalues, In the same manner probabilities for λ y = ±1/ √ 2 are 0.5, while λ z = 1/ √ 2 is fixed (probability 1). Our GDTWA prescription thus reduces to the DTWA sampling introduced in Ref. [51].

III. COMPARISON WITH EXACT DIAGONALIZATION
We now proceed to benchmark the validity of the GDTWA by comparison with the dynamics obtained from numerical exact diagonalization (ED). Specifically, we will consider the dynamics induced by dipolar interactions, which are currently investigated in a variety of experimental platforms, including polar molecules [52], magnetic atoms [29,34], NV centers [53], etc. The Hamiltonian for such systems can be written in the form of an anisotropic XXZ spin model: where V i,j = V dd (1 − 3 cos 2 θ i,j )/r 3 i,j is the strength of the interaction between two atoms i and j, whose distance is r i,j and θ i,j is the angle between the vector r i,j and the dipole orientation typically set by an external quantization field. For magnetic atoms such as chromium (S = 3) and erbium (S = 19/2 for fermions, and S = 6 for bosons), V dd = µ 0 (µ B g) 2 /4π, with µ 0 the magnetic permeability of vacuum, µ B the Bohr magneton, and g the Landé g-factor. For convenience, we define a bare dipolar interaction strength V = V dd /d 3 , with d the smallest spacing between different lattice sites.
Spin-level Populations: An appealing property of magnetic atoms is their sizable spin, which by itself is responsible for enhanced dipolar interactions. A large spin provides a great degree of tunability through the control of the various multiple internal spin states. In recent experiments the spin dynamics has been activated by preparing two types of initial states: i) a spin coherent state of atoms [17,29] |ψ(t = 0) = ⊗ j exp(iθŜ j y ) |−S j , with θ ∈ [0, π] the tipping angle, or ii) all atoms in the same Zeeman level [31][32][33][34]   where m 0 = −S, −S + 1, . . . , S, and j = 1, 2, . . . , N labels atoms (see also A for discussions on initial states in GDTWA). In Fig. 2, we compute the evolution of population in different Zeeman levels m averaged over all spins, n m (t) = (1/N ) j | ψ(t)|m j | 2 , starting from the different initial states for both scenarios with S = 3 and S = 19/2. We compare a simulation using the GDTWA with the exact ED result. The very good agreement between GDTWA and ED shows that the GDTWA captures the dynamics for all Zeeman levels well under dipolar exchange. While, as expected the GDTWA provides quantitatively exact results for short times, remarkably it features also qualitatively the long-time behavior well. Note that an attempt to simulate the population dynamics of Fig. 2 with a TWA approach with Gaussian sampling of three spin-variables leads to unphysical results as shown in B.

Single-spin Density Matrix:
In addition to Zeeman level populations and spin projections Ŝ i x,y,z , the GDTWA also provides access to all other elements of any single-spin density matrix element. It is important to note that this cannot be obtained with a conventional Gaussian sampling. In particular, while in principle a state-tomography can also be performed from combinations of expectation values of powers of spin operators, the Gaussian sampling does not properly account for high order moments. For example, in the case of a spin coherent state along the z-axis the Gaussian sampling provides incorrect results for expectation values Ŝ n x,y with n > 3. This implies that for models with S > 3/2, a state tomography leads to significant errors for density matrix elements, already in the initial state (see B). In the GDTWA, the ability to predict all density matrix elements allows us to compute entanglement measures such as the Renyi entropy for a single spin in the coupled system, S c α = 1 1−α log 2 Tr[ρ α c ], whereρ c is the reduced density matrix of a single spin c, and α ≥ 0. In Fig. 3, we com-  ) and GDTWA (dots). The system is a 1D chain with N = 8 atoms with spin S = 2 and with dipolar interactions (dipole orientation perpendicular to the chain). The initial state is a spin coherent state polarized along the chain direction, i.e. of type i) with θ = π/2. The GDTWA captures the interaction induced spinsqueezing (ξR < 1) nearly perfectly.
pare the evolution of the second Renyi entropy (α = 2) computed with GDTWA and ED for the same scenarios as in Fig. 2. For both initial states, with S = 3 and S = 19/2 respectively, we find excellent agreement for all timescales. Note that in C, we provide further comparisons for type i) initial state with different θ to better analyze the validity of the GDTWA. There we observe that the GDTWA can provide nearly exact predictions, especially in cases where the system quickly reaches a steady state with high local entropy. Nevertheless, deviations from the exact Renyi entropy evolution are seen at low angles (θ → 0). This is an interesting observation that needs further investigation given that the θ → 0 limit is where a simple mean-field treatment is enough to accurately reproduce the short time dynamics.
Spin Squeezing: In a many-body system, interactions can generate quantum spin squeezing, which can fundamentally improve measurement precision [54][55][56]. Spin squeezing is of interest to many current experiments [25,[57][58][59]. Here we demonstrate that the GDTWA is also applicable for studying such a phenomenon. In particular, we study the spin squeezing parameter first introduced by Wineland et al. as a measure of phase sensitivity [60,61], where ∆S 2 ⊥ = Ŝ 2 ⊥ − Ŝ ⊥ 2 is the variance measured perpendicular to the collective spin S direction, and | S| = Ŝx 2 + Ŝy 2 + Ŝz 2 is the length of the collective spin. As indicated by this definition, the spin squeezing parameter involves quantum correlations among many atoms. For a spin coherent state, ∆S 2 ⊥ = N S/2, ξ 2 R = 1.
On the other hand ξ 2 R < 1 implies reduced phase noise. As shown in Fig. 4, the result obtained with GDTWA again agrees almost perfectly with those obtained from ED. Note that compared to Eq. (18), here we use a slightly modified method to compute collective spinobservables such as (Ŝ x ) 2 . In various benchmark situations we find that this modified procedure provides much more precise predictions for the spin-squeezing than Eq. (18). We discuss this in D.

IV. SPIN DYNAMICS IN MANY-BODY SYSTEMS
While the system sizes accessible with ED are remarkably small for large S, the GDTWA enables us to overcome this limitation and to perform investigations of spin dynamics in large ensembles with large S. This capability is very important to quantitatively perform comparisons to experiments, in particular in the presence of collective behaviors, where finite size effects matter. In the following sections we will employ the GDTWA to obtain new insight for those systems by studying: time-dependent re-distribution of populations between Zeeman levels and resulting quantum thermalization for this observable, the evolution of entanglement among the dipoles, the evolution of the contrast of spin coherence, and spin-squeezing.

A. Spreading of Zeeman level populations
Here, we apply the GDTWA approach to study the spin dynamics in a 3D lattice of dipoles for various S, with the interaction Hamiltonian given by Eq. (21). We prepare the system in the type ii) initial state, which is not an eigenstate of the many-body Hamiltonian H dd , and thus the population of the initial Zeeman level will generally change with time. Specifically, the dipolar exchange will lead to a redistribution of population over all 2S + 1 levels (see Fig. 5(a)). As demonstrated in Ref. [34], the effective dipolar exchange strength given by V γ(S, m 0 ) can be used as the characteristic energy scale of the spin dynamics, with γ(S, m 0 ) given by where m 0 is the initially populated spin level. Correspondingly we defineṼ = V γ(S, m 0 ). As shown in Fig. 5(b), for various values of S, the initial dynamics happens at a characteristic rate 1/Ṽ .
It is worth to note that for these initial conditions within a mean-field approximation one observes the spinlevel populations to remain fixed at the initial value without any evolution. Quantum fluctuations of the initial state are what drives the dynamics in this case and therefore need to be accounted for. The GDTWA includes those initial fluctuations in an exact way and thus can describe the evolution of the level populations. Note that these simulations are done for a 3D system with ∼ 200 atoms, which corresponds to a full Hilbert space ∼ 10 240 states for S = 6 and is far beyond the capacity of ED.
At long times, when the system thermalizes, the spin dynamics shows saturation to different values for different S. According to the theory of closed system thermalization, given e.g. by the Eigenstates Thermalization Hypothesis (ETH), it is expected that for a nonintegrable system without disorder [62][63][64][65], at long times the expectation values of local observables can be effectively described by a thermal equilibrium state,ρ ∞ = exp(−βĤ T − µŜ z )/Ξ, with Ξ = Tr[exp(−βĤ T − µŜ z )], H T =Ĥ dd , andŜ z = iŜ i z . The effective inverse temperature β and chemical potential µ are set by the total energy and magnetization, respectively, which are the conserved quantities in our problem. Those are determined by the initial state and conserved under the unitary evolution. In our system, the total energy for the initial state is E 0 = i,j =i V i,j m 2 0 /2. In a high temperature limit β → 0, energy conservation leads to For a 3D lattice where dipolar interactions vary in sign, E 0 ≈ 0 and thus we take β = 0 in the following. Then µ can be easily determined from As shown in Fig. 5 (c) and (d), the long-time populations obtained from the GDTWA agree extremely well with the quantum thermalization predictions for both integer and half-integer S. For integer S, an initial state with m 0 = 0 results in equal population on all Zeeman levels in equilibrium. This is not the case for half-integer spins where the steady state populations are not equal. The different steady-state behavior seen in integer and halfinteger spin systems is reminiscent of their different low energy spectra at equilibrium. The agreement shown here suggests that in the longtime limit GDTWA simulations for simple single-spin observables approach the results expected from closed system quantum thermalization, an effect also observed previously for different initial states for S = 1/2 [48] and S = 3 [29] models. Note that this implies that in models which feature quantum thermalization, i.e. in models which are not many-body localized, and with local observables that thermalize quickly, the GDTWA can give accurate predictions for simple observables not only at short times, but over all relevant timescales.
Furthermore, the knowledge of the equilibrium values also allows us to use the GDTWA to investigate how such a quantum system approaches thermalization. Recent works have suggested a power-law relaxation in the  long-time dynamics [66,67]. Interestingly, as illustrated in the inset of Fig. 5 (b), here we observe a power-law decay towards thermalization at long times, which roughly follows n m0 (t) − n m0 (∞) ∼ 1/t 2 , with n m0 (t) the population on state m 0 at time t.

B. Entanglement build-up measured through the Renyi Entropy
As mentioned in Sec. III, in GDTWA simulations we can easily compute the second Renyi entropy S c 2 = − log 2 Tr[ρ 2 c ] for a subsystem of a single spin. The entropy is quantifying the purity of the subsystem and thus provides useful information on the bi-partite entanglement with all other spins, given that the state of the full system remains a pure state. In Fig. 6 we show the second Renyi entropy for a single spin during the spin dynamics of Fig. 5(b), withρ c chosen as the reduced density matrix of a central spin in the cube. Since our initial state is a pure state with S c 2 = 0, an increasing S c 2 (t) > 0 indicates a build-up of entanglement in the system. For all values of S, S c 2 increases quickly at short-times in an algebraic way. It gradually approaches the maximum second Renyi entropy of log 2 (2S + 1) which is allowed by the local Hilbert space size of 2S +1, and which increases with S. This suggests that, atoms with large S, rather than behaving more classically, feature richer quantum effects due to inter-spin entanglement. And since for dipolar systems, such as the magnetic atoms discussed above, larger S is associated with a larger net interaction strength, this also implies a faster generation of entanglement in the spin dynamics.

C. Dynamics of contrast and intra-spin correlations
Another useful probe of many-body effects is the contrast of spin coherence, defined as C S = Ŝx 2 + Ŝy 2 /(N S), withŜ x,y,z = jŜ x,y,z j . It has been measured via Ramsey spectroscopy in many experiments [15,68,69]. Fig. 7(a) shows the evolution of C S under the dipolar interactions in Hamiltonian (21), starting from the type i) spin coherent initial state with tipping angle θ = π/2. Using a short-time expansion analysis, we find that the contrast initially decays as [70] Compared with the dynamics discussed in the previous section, where a larger spin S can provide an enhanced dipole-dipole interaction ∝Ṽ ∝ S 2 V , here the contrast decays slower than the timescale ∝ 1/Ṽ . Instead, the characteristic interaction strength for the contrast dynamics isṼ c = √ SV . This is seen in Fig. 7(a), where as a function ofṼ c t all lines collapse onto a single one for short times.
As mentioned in Sec. II, such an initial spin state can also be represented in a conventional Gaussian TWA for three spin-variables. However, for S > 1/2, quantum intra-spin correlations can develop [71]. As discussed in Sec. III, since in a Gaussian TWA only the spin op-eratorsŜ j x,y,z are involved, such intra-spin correlations generally cannot be well captured. To illustrate this effect, in Fig. 7(b) we compare the evolution of j (Ŝ j z ) n (n = 4, 6), computed with GDTWA and Gaussian TWA for three spin-variables. The Gaussian TWA not only provides incorrect initial values, but also shows an increased deviation as the spin dynamics proceeds. To reduce finite-size effect, we consider a 1D chain with various large N to adjust for a fixed N S = 240. Dipoles are aligned perpendicular to the chain axis (θi,j ≡ π/2). The initial state is a spin coherent state polarized along the chain direction, i.e. of type i) with θ = π/2.

D. Quantum spin squeezing
Lastly, to further characterize the quantum correlations generated during the dipolar spin dynamics, we study the spin squeezing parameter ξ 2 R introduced in Sec. III for ensembles of spins with various S. It is straightforward to show that ξ 2 R ≥ 1/(2N S). Hence, to make a fair comparison for systems of different spin S, in Fig. 8, we keep N S constant and calculate ξ 2 R for various S. While for S > 1/2 spin squeezing is no longer an entanglement witness given that part of the contributions come from intra-spin correlations [71], nevertheless the calculations do suggest an enhanced phase sensitivity potentially provided by large S particles. Explicitly, compared to the conventional spin-1/2 particles, a system of S = 2 shows significantly larger spin squeezing and also the timescale for reaching maximum squeezing is shortened, which could be a favorable feature in the presence of dephasing.

V. CONCLUSION AND OUTLOOK
The numerical approach discussed here is capable of capturing beyond mean-field effects and entanglement build-up efficiently. For the case of dipolar interactions presented here we have shown that the GDTWA gives an excellent approximation for short times and features correct quantum-thermalization in the long-time limit. This implies that in models which thermalize, e.g. which are not many-body localized because of disorder, the method can give accurate predictions over all time-scales. The method can be applied to any model consisting of discrete coupled Hilbert-spaces, such as spin systems featuring SU (N ) super-exchange interactions [72], the so called Subir-Ye-Kitaev models [73][74][75], Bose-Hubbard models [76], etc. It furthermore allows one to systematically treat the case where one builds clusters of spin-1/2 particles in order to better account for quantum correlations. By comparing the dynamics obtained for different cluster sizes, one could benchmark the convergence of the phase space method in regimes where no other methods are available for verification. This approach was recently suggested in Ref. [42], however without utilizing initial discrete distributions.
For the time-evolution in our GDTWA approach, correlations of GGMs such as λ i µ (0)λ j ν (0) , λ i µ (0)λ j ν (0)λ k κ (0) are important. Here, we will show that for diagonal states the quantum mechanical correlations of the initial state are perfectly reproduced by our discrete distributions from Eq. (15). Correlations between sites trivially factorize since we assume a product state, and the lack of inter-site correlations is correctly re-produced by our discrete distributions (15), which also factorize between sites. Here, we thus focus on a single spin (site indices are dropped for simplicity). By diagonal state |α 0 with α 0 = 1, 2, . . . , D, we mean any eigenstate of the diagonal GGMs defined in Eqs. (10), which are eigenstates ofŜ z .
We first note that for diagonal states, any power of a single GGM is fully re-produced by an average over our distribution since forρ = |α 0 α 0 |, where we used the expansion of the GGM into the orthonormal eigenbasis Λ µ = aµ a µ |a µ a µ |. Since the discrete distributions defined in Eq. (15) are furthermore independent for different µ, correlations between different GGMs always factorize: with arbitrary powers n 1 , n 2 , . . . , n D ≥ 0. We will show that this properties also holds for correlations for the diagonal quantum state. We proceed by showing that for µ = ν The symmetrization super-operator, S, is necessary since generally different GGMs cannot be measured simultaneously as they do not always commute. Note that (A3) is always true if |α 0 is an eigenstate of Λ µ and Λ ν . This is trivially the case for all diagonal matrices, and therefore for all diagonal GGMs from Eqs. (10). Moreover, note that for an arbitrary diagonal state |α 0 there are only 2(N − 1) GGMs for which |α 0 is not an eigenstate, with ξ = 1, 2, . . . , α 0 − 1, α 0 + 1, . . . , N . For all other off-diagonal GGMs, |α 0 is an eigenstate with zero eigenvalue. We therefore have to only show that Eq. (A3) also holds if at least one of the matrices Λ µ and Λ ν is from Eqs. (A4). It is straightforward to compute that for µ = ν Re α 0 | R α0 ξ I α0 ξ |α 0 = Re whereΛ µ is any GGM not from (A4). Therefore, (A3) is correct. Additionally, we can use the fact that any power of a GGM from Eqs. (8), Eqs. (9), and Eqs. (10) is either a diagonal matrix and trivially has |α 0 as eigenstate, or it is proportional to itself, therefore also (µ = ν) holds for arbitrary integer powers of n µ and n ν . Lastly, it's straightforward to see that the pairwise property (A9) also leads to a factorization of any other symmetrized product of multiple GGMs. To realize this, consider that Here we used the group property of the GGMs, i.e. that the product of any number of GGMs is again proportional to some other GGM, Λ ξ . If Λ ξ commutes with (Λ κ ) nκ we do not need to symmetrize, otherwise κ = ξ and we can again use Eq. (A9) to see that α 0 | SΛ ξ (Λ κ ) nκ |α 0 factorizes. From a repeated application of this argument it follows that In particular, for states |ψ that are not diagonal, but can be obtained from a diagonal state via any singlesite unitary transformation (Û †Û = 1), |ψ =Û |α , we can transform the Hamiltonian into the rotated basiŝ H =Û †ĤÛ , and then obtain the equations of motion for GGMs in a similar way as Eq. (14), withĤ replaced byĤ . The dynamics of an observable of interestÔ can the be found from evaluatingÔ =Û †ÔÛ . Since in this rotated basis the initial state becomes diagonal, the expectation values ofÔ can be correctly described with the discrete sampling, and thus ψ|Ô |ψ in the original basis can be obtained. For the type i) initial states (used in the main text) with arbitrary θ, we can useÛ = e iθŜ y and perform sampling in the rotated basis with initial state |−S . For initial states that are a coherent superposition of several Zeeman levels, |ψ = m c m |m , with m = −S, −S + 1, ...S, c m real numbers and m c 2 m = 1, the unitary transformation can be constructed from m ]Â −S,m , witĥ A −S,m = |−S m| + |m −S| hermitian operators. The type ii) initial states considered in the main text constitute a special case in this category. In the case that initial states vary among particles, the overall transformation can be straightforwardly constructed accounting for spatial inhomogeneities:Û N = ⊗ jÛj withÛ † jÛ j = 1 j .  In Fig. 3(a) we have shown that for S = 3 and tipping angle θ = π/2 the GDTWA approach captures the second Renyi entropy nearly perfectly. Here to further show the validity of our approach, we provide results for several other tipping angles in Fig. 11. Except for θ = 0.2π, these comparisons show a good agreement between the GDTWA and ED, suggesting that generally for θ not too small, the GDTWA approach works well for calculating the second Renyi entropy. Surprisingly, in the small θ case, in which even a pure mean-field simulations becomes valid to describe the short time dynamics, due to generally small build-up of entanglement, the GDTWA shows some intrinsic production of entanglement which is unphysical. We attribute this to some incorrect correlation build-up in higher-order correlations, which are incorrectly captured due to the TWA being a lowest order expansion (in ). Indeed this is corroborated by Fig. 11(b), where we show the case of θ = 0. Quantum mechanically, one would not expect any dynamics in this case, since the initial state is an eigenstate of the Hamiltonian. However, the GDTWA still shows some dynamics. The deviation (difference with zero Renyi entropy) grows initially ∝ (V t) 4 , which indicates that higher order corrections to the GDTWA are relevant. It is again worth pointing out, that nevertheless in the case of large entanglement build-up, these corrections seem to contribute only a relatively minor correction. This means that specifically in the regime where mean-field simulations fail, the GDTWA can provide very good estimations.
Appendix D: Calculating spin squeezing in the GDTWA To calculate the spin squeezing parameter ξ 2 R , we need to evaluate the expectation value of linear and quadratic operators such asŜ x and (Ŝ x ) 2 . In the GDTWA, we first evaluate the expansion coefficients for the general operators in the GGM basis [Eq. (17)]: Then, we approximate the quantum expectation values from averages over the classical trajectories: [Eq. (18)]: Alternatively, in the expansion of (Ŝ x ) 2 in Eq. (D2) one could also ignore the distinct treatment given to the diagonal terms, and instead compute the expectation value (Ŝ x ) 2 as In principle the first term in Eq. (D4) takes the singlespin part of the collective observable into account more precisely. Note, e.g. that for S = 1/2 it correctly reproduces the constant value for σ 2 x = 1 = 1. Nevertheless, by comparisons with ED calculations, we find that in practice, especially for finite range interactions, Eq. (D5) describes the evolution of ξ 2 R quantitatively better than Eq. (D4). This is illustrated in Fig. 12, where a comparison with ED shows a nearly perfect agreement using Eq. (D5), while a clear deviation on the logarithmic scale is visible when using Eq. (D4). Therefore, in the main text we adopted Eq. (D5) for evaluating spin squeezing.