# A generalized phase space approach for solving quantum spin dynamics

###### Abstract

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 . We show that the GDTWA can accurately simulate dynamics of large ensembles in arbitrary dimensions. We apply it for 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 can correctly reproduce 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 models which feature quantum-thermalization at long times. We apply our method to study dynamics in large spin-models and compute experimentally accessible observables such as Zeeman level populations, contrast of spin coherence, spin squeezing, and entanglement quantified by single-spin Renyi entropies. We reveal that large systems can feature larger entanglement than corresponding 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 out-of-equilibrium dynamics of large isolated quantum many-body models in controlled environments Bloch (2018); Cirac and Zoller (2012); Bloch et al. (2012); Blatt and Roos (2012). 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 Orús (2014); Schollwöck (2011); Verstraete et al. (2008); Vidal (2003); White and Feiguin (2004); Daley et al. (2004)), particular ansatz wave-functions (e.g. in combination with variational Monte-Carlo evolution Carleo et al. (2012); Cevolani et al. (2015); Ido et al. (2015); Carleo and Troyer (2017)), or dilute systems (e.g. by making use of a clusterization Hazzard et al. (2014a); Piñeiro Orioli et al. (2018)).

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 have been implemented with many different setups, e.g. using polar molecules Yan et al. (2013); Hazzard et al. (2014a), Rydberg atoms Labuhn et al. (2016); Zeiher et al. (2017); Bernien et al. (2017); Barredo et al. (2018); Guardado-Sanchez et al. (2018); Takei et al. (2016); Piñeiro Orioli et al. (2018), trapped ions Gärttner et al. (2017); Bohnet et al. (2016); Neyenhuis et al. (2017) and cavity QED systems Norcia et al. (2018); Davis et al. (2019), recently also models with larger spins have become a research focus in particular for experiments with magnetic atoms Lepoutre et al. (2019); Baier et al. (2016); de Paz et al. (2016a, 2013, b); Patscheider et al. (2019). The large spin degrees of freedom in systems poses a much more stringent requirement for numerical treatment compared to systems. The full Hilbert space involves states for particles, which renders exact diagonalization (ED) impossible already for small . 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 Zaletel et al. (2015); Haegeman et al. (2016); Schachenmayer et al. (2013). 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 . The method is based on the well known truncated Wigner approximation (TWA) Graham (1973); Steel et al. (1998); Blakie et al. (2008); Polkovnikov (2010) 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 spin-variables, but all density matrix elements for each spin Lepoutre et al. (2019); Wurtz et al. (2018); Patscheider et al. (2019); Davidson and Polkovnikov (2015). 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 our generelized discrete TWA approach (GDTWA), reduces to the previously proposed discrete TWA method (DTWA) Schachenmayer et al. (2015a), which has been remarkably succesful in predicting model dynamics Schachenmayer et al. (2015b); Pucci et al. (2016); Piñeiro Orioli et al. (2017); Acevedo et al. (2017); Piñeiro Orioli et al. (2018); Czischek et al. (2018).

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 Lepoutre et al. (2019); de Paz et al. (2016a, 2013, b) and Er atoms Patscheider et al. (2019)) 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

### ii.1 Phase space sampling for quantum spin systems

The Hamiltonian for a system consisting of spin- particles coupling to each other via two-body interactions can be written as

(1) |

with , and . 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 particle, , can be obtained via the Heisenberg equations of motion (we set throughout this paper):

(2) |

These equations generally depend on high order inter-spin correlations such as , , , for and in practice it is impossible to solve them exactly for large . In a mean-field ansatz one assumes that such spin correlations factorize and can be written in terms of single particle observables, , e.g., . Then Eqs. (2) turn into coupled closed equations which can be easily solved numerically, and which correspond to equations for classical spin-variables^{1}^{1}1Note that more general models than Eq. (1) for can also depend on intra-spin terms such as quadratic fields . 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.. 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 hierarchies Pucci et al. (2016). However, then already retaining second-order correlations increases the complexity of Eqs. (2) to and furthermore higher order corrections can typically lead to numerical instabilities. An alternative approach to simulate quantum correlations while retaining a complexity is to resort to a phase space description of the quantum system Blakie et al. (2008); Polkovnikov (2010) 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 and , respectively, and the Weyl symbol is denoted as . The expectation value of any operator can then be computed as an integral in phase space:

(3) |

Here, is the Wigner function corresponding to the Weyl symbol of the density matrix : , 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:

(4) |

Here, and are classical trajectories of and , respectively, with and . 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 Polkovnikov (2010); Blakie et al. (2008). 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 corresponds to the observable Polkovnikov (2010).

For a system of coupled spin- particles, a natural way to formulate the TWA dynamics consists of replacing , and using the spin-variables (with and ), corresponding to the Weyl symbols of , as phase space variables. The trajectories correspond to those obtained from the mean-field approximation of Eq. (2) with . For particular initial states, the Wigner function can be easily found. Taking for example a product state , where each spin is described by a spin coherent state, then for large the Wigner function can be well approximated by a Gaussian. In the case of a state initially polarized along the direction, the Wigner function factorizes, and for each spin takes the simple form

(5) |

where denotes the delta function. For this state the phase space value of is determined as , while the values of () fluctuate according to a Gaussian distribution with zero mean and a variance (see Fig. 1). This width of the distribution of the classical variables reflects the uncertainty relation intrinsic to the quantum mechanical operators and . For large , , suggesting vanishing effects from quantum fluctuations.

For such initial states, the TWA evolution now reads

(6) |

where the Weyl symbol corresponds again to symmetrized observables. For example, for an inter-spin correlation such as at distinct sites , simply . 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 spin-component for a single spin can now depend on the initial conditions of all other spin-variables via some function, e.g. . Therefore, for observables computed in the TWA, such as with ,

(7) |

and correlations can build up with time. Generally, the factorization in the last equation would only hold if would only depend on the variables of the spin itself.

### ii.2 Generalized discrete Truncated Wigner approximation (GDTWA)

There are some important shortcomings of the Gaussian TWA approach with spin-variables introduced in the previous section. Most importantly, the scheme is tailored to problems relying on the spin operators, both in the Hamiltonian and for measurements. Most experimental scenarios go beyond this limitation, by including e.g. quadratic fields and measurements of Zeeman state populations Lepoutre et al. (2019); Patscheider et al. (2019). Furthermore, the approach is insufficient for a full state-tomography in large 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 spin-variables to elements of higher-dimensional generalizations of Bloch vectors Bertlmann and Krammer (2008). This can be accomplished by noticing that for a spin- atom with spin states, its density matrix consists of elements. Correspondingly, we can define hermitian operators, (hats are dropped for clarity of notation), using the generalized Gell-Mann matrices (GGM) and the identity matrix ^{2}^{2}2Note that for convenience, we use a slightly different normalization factor than in the definition of usual (generalized) Gell-Mann matrices Bertlmann and Krammer (2008), which will simplify operator expansions.:

(8) | ||||

(9) | ||||

(10) | ||||

(11) |

For each spin , these matrices are orthonormal, , and constitute a complete local basis. Hence any single-spin operator can be represented via

(12) |

and . Consequently, more general Hamiltonians than Eq. (1), with one- and two-body terms can now be represented as

(13) |

with and . The evolution of expectation values of GGMs can now be computed as

(14) |

For each GGM we define a corresponding real phase space variable, . 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 spin-variables discussed after Eq. (2), our classical equations of motions follow from assuming a factorization between different sites for any nonequal site-index in Eq. (14), and by replacing .

To define a phase space probability distribution for the inital state, we decompose each via its eigenvectors with corresponding eigenvalues , . As in spin-1/2 systems, where the Pauli matrices can be measured to be in a projective measurement, the correspond to possible measurement outcomes of . We will focus on the experimentally relevant case of initial product states . Then, for each we can define a corresponding probability distribution for our phase space variables, which we limit to a discrete set of values

(15) |

The overall distribution factorizes for different variables on the same site and between sites, such that the probability for a configuration of all being a certain combination of the eigenvalues for and is given by .

From the distribution (15), we can compute the expectation value of any observable on a fixed site via the expansion (12),

(16) |

Here we defined the notation , denoting the statistical average over any combination of eigenvalues for the 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 spin-operators on a site , and do not form a complete local basis for any operator. Thus, for example the observable cannot be expanded as linear superposition of the matrices and has to be evaluated as 4-th order moment from the probability distribution. As a consequence, for a spin coherent state polarized along , the Gaussian distribution does not reproduce 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 for all and . 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 . Here is a vector of all phase-space variables , and the matrix depends in general on the configuration at time . We solve those equations numerically for the selected initial condition and obtain a mean-field trajectory, . 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 at time as .

If the problem is linear, i.e. if is time-independent, e.g. if there are only single-spin terms in Eq. (14), then the time-dependent solution in our scheme is given by , and since via Eq. (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 will also depend on intra-spin correlations in the initial state, also within a single spin, e.g. , 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. with . This is not necessarily true for the quantum mechanical intra-spin correlations of a generic single-spin state. For example, for some states , , 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 as where , and 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 , 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 Wurtz et al. (2018); Davidson and Polkovnikov (2015), i.e. for all 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 Schachenmayer et al. (2015c).

To compute the time-dependent expectation value of an arbitrary multi-spin observable, , we expand

(17) |

and approximate

(18) |

This equation is a discrete analog of the TWA method from Eq. (6) for a phase space extended to a high-dimensional generalized Bloch sphere. Approximating Eq. (18) with a finite number of samples, , 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, . After selecting a single random configuration for all phase space variables according to , we denote the subsequent classical evolution of the variables at site for this sample “” as . Then, the evolution of is computed as

(19) |

Note that in contrast to a continuous Wigner function approach, here in principle only a finite number of numerical samples is needed to compute Eq. (18). In practice, however, this number increases exponentially with the . The 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 samples sufficient. Note that we find that for collective observables the number of samples necessary for convergence typically decreases with .

To further explain the main idea behind this approach and its capabilities, we first consider as an example an array of spin particles. In this case, the three nontrivial for each spin are proportional to standard Pauli matrices, , each of which has eigenvalues , and eigenstates described by states aligned (anti-aligned) along the corresponding direction: . For a spin initially polarized as , the initial values of in phase space take its two eigenvalues, , with equal probability . In the same manner probabilities for are , while is fixed (probability 1). Our GDTWA prescription thus reduces to the DTWA sampling introduced in Ref. Schachenmayer et al. (2015c).

For , there are three states for each spin, and consists of the identity matrix and 8 nontrivial Gell-Mann matrices:

(20) |

An initially polarized state , can now be represented via probabilities for eigenvalues of . Specifically, with probability one chooses at value , and with probability one chooses the value . In the following, we denote this as , with probability . Then, further we have fixed with , with , and with (See Fig. 1 for an illustration of this distribution).

While the initial states discussed above are still spin coherent states, the state , however, is an example of a state that cannot be represented with the conventional Gaussian TWA approach. In the GDTWA it can be easily simulated with , , , , , , and , .

## 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 Bohn et al. (2017), magnetic atoms Lepoutre et al. (2019); Patscheider et al. (2019), NV centers Choi et al. (2017), etc. The Hamiltonian for such systems can be written in the form of an anisotropic XXZ spin model:

(21) |

where is the strength of the interaction between two atoms and , whose distance is and is the angle between the vector and the dipole orientation typically set by an external quantization field. For magnetic atoms such as chromium () and erbium ( for fermions, and for bosons), , with the magnetic permeability of vacuum, the Bohr magneton, and the Landé g-factor. For convenience, we define a bare dipolar interaction strength , with 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 Yan et al. (2013); Lepoutre et al. (2019) , with the tipping angle, or ii) all atoms in the same Zeeman level de Paz et al. (2016a, 2013, b); Patscheider et al. (2019) , where and 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 averaged over all spins, , starting from the different initial states for both scenarios with and . 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 , 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 -axis the Gaussian sampling provides incorrect results for expectation values with . This implies that for models with , 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, , where is the reduced density matrix of a single spin , and . In Fig. 3, we compare the evolution of the second Renyi entropy () computed with GDTWA and ED for the same scenarios as in Fig. 2. For both initial states, with and 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 (). This is an interesting observation that needs further investigation given that the 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 Kitagawa and Ueda (1993); Ma et al. (2011); Pezzè et al. (2018). Spin squeezing is of interest to many current experiments Esteve et al. (2008); Leroux et al. (2010); Bohnet et al. (2016); Hosten et al. (2016). 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 Wineland et al. (1992, 1994),

(22) |

where is the variance measured perpendicular to the collective spin direction, and 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, , . On the other hand 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 spin-observables such as . 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 , the GDTWA enables us to overcome this limitation and to perform investigations of spin dynamics in large ensembles with large . 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.

### iv.1 Spreading of Zeeman level populations

Here, we apply the GDTWA approach to study the spin dynamics in a 3D lattice of dipoles for various , 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 , 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 levels (see Fig. 5(a)). As demonstrated in Ref. Patscheider et al. (2019), the effective dipolar exchange strength given by can be used as the characteristic energy scale of the spin dynamics, with given by

(23) | ||||

where is the initially populated spin level. Correspondingly we define . As shown in Fig. 5(b), for various values of , the initial dynamics happens at a characteristic rate .

It is worth to note that for these initial conditions within a mean-field approximation one observes the spin-level 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 atoms, which corresponds to a full Hilbert space states for 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 . According to the theory of closed system thermalization, given e.g. by the Eigenstates Thermalization Hypothesis (ETH), it is expected that for a non-integrable system without disorder D’Alessio et al. (2016); Deutsch (2018); Basko et al. (2006); Nandkishore and Huse (2015), at long times the expectation values of local observables can be effectively described by a thermal equilibrium state, , with , , and . 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 . In a high temperature limit , energy conservation leads to

(24) |

For a 3D lattice where dipolar interactions vary in sign, and thus we take in the following. Then can be easily determined from , with , and we can obtain the equilibrium population in a spin state:

(25) |

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 . For integer , an initial state with 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 half-integer spin systems is reminiscent of their different low energy spectra at equilibrium.

The agreement shown here suggests that in the long-time 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 Acevedo et al. (2017) and Lepoutre et al. (2019) 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 Santos and Torres-Herrera (2017); Borgonovi et al. (2019). 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 , with the population on state at time .

### iv.2 Entanglement build-up measured through the Renyi Entropy

As mentioned in Sec. III, in GDTWA simulations we can easily compute the second Renyi entropy 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 chosen as the reduced density matrix of a central spin in the cube. Since our initial state is a pure state with , an increasing indicates a build-up of entanglement in the system. For all values of , increases quickly at short-times in an algebraic way. It gradually approaches the maximum second Renyi entropy of which is allowed by the local Hilbert space size of , and which increases with . This suggests that, atoms with large , 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 is associated with a larger net interaction strength, this also implies a faster generation of entanglement in the spin dynamics.

### iv.3 Dynamics of contrast and intra-spin correlations

Another useful probe of many-body effects is the contrast of spin coherence, defined as , with . It has been measured via Ramsey spectroscopy in many experiments Hazzard et al. (2014a); Martin et al. (2013); Bromley et al. (2018). Fig. 7(a) shows the evolution of under the dipolar interactions in Hamiltonian (21), starting from the type i) spin coherent initial state with tipping angle . Using a short-time expansion analysis, we find that the contrast initially decays as Hazzard et al. (2014b)

(26) |

Compared with the dynamics discussed in the previous section, where a larger spin can provide an enhanced dipole-dipole interaction , here the contrast decays slower than the timescale . Instead, the characteristic interaction strength for the contrast dynamics is . This is seen in Fig. 7(a), where as a function of 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 , quantum intra-spin correlations can develop Vitagliano et al. (2014). As discussed in Sec. III, since in a Gaussian TWA only the spin operators are involved, such intra-spin correlations generally cannot be well captured. To illustrate this effect, in Fig. 7(b) we compare the evolution of (), 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.

### iv.4 Quantum spin squeezing

Lastly, to further characterize the quantum correlations generated during the dipolar spin dynamics, we study the spin squeezing parameter introduced in Sec. III for ensembles of spins with various . It is straightforward to show that . Hence, to make a fair comparison for systems of different spin , in Fig. 8, we keep constant and calculate for various . While for spin squeezing is no longer an entanglement witness given that part of the contributions come from intra-spin correlations Vitagliano et al. (2014), nevertheless the calculations do suggest an enhanced phase sensitivity potentially provided by large particles. Explicitly, compared to the conventional spin-1/2 particles, a system of 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 super-exchange interactions Gorshkov et al. (2010), the so called Subir-Ye-Kitaev models Sachdev and Ye (1993); Maldacena and Stanford (2016); Schmitt et al. (2019), Bose-Hubbard models Nagao et al. (2019), etc. It furthermore allows one to systematically treat the case where one builds clusters of spin- 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. Wurtz et al. (2018), however without utilizing initial discrete distributions.

### Acknowledgements

We thank Francesca Ferlaino’s Er group and Bruno Laburthe-Tolra’s Cr group and Chunlei Qu for fruitful discussions. We also thank Asier Pineiro and Itamar Kimchi for reviewing the manuscript. Funding: B.Z. is supported by the NSF through a grant to ITAMP. A.M.R is supported by the AFOSR grant FA9550-18-1-0319 and its MURI Initiative, by the DARPA and ARO grant W911NF-16-1-0576, the DARPA DRINQs grant, the ARO single investigator award W911NF-19-1-0210, the NSF PHY1820885, NSF JILA-PFC PHY1734006 grants, and by NIST. J.S. is supported by the French National Research Agency (ANR) through the Programme d’Investissement d’Avenir under contract ANR-11-LABX-0058_NIE within the Investissement dâAvenir program ANR-10-IDEX-0002-02, by IDEX project “STEMQuS”, and through computing time at the Centre de calcul de l’Université de Strasbourg.

## Appendix A Initial states that can be represented within the GDTWA

For the time-evolution in our GDTWA approach, correlations of GGMs such as , 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 with , we mean any eigenstate of the diagonal GGMs defined in Eqs. (10), which are eigenstates of .

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 ,