A generalized phase space approach for solving quantum spin dynamics

TL;DR

The paper introduces GDTWA, a new semi-classical phase-space sampling method for simulating quantum spin dynamics in large, interacting lattice systems with arbitrary spin, capturing complex thermalization and entanglement phenomena.

Contribution

It presents GDTWA, a novel numerical approach that accurately models quantum spin dynamics beyond mean-field approximations for large, high-spin systems with long-range interactions.

Findings

GDTWA accurately reproduces short- and long-time dynamics.

The method captures quantum-thermalization effects.

Large S systems can exhibit greater entanglement than S=1/2 systems.

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 $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.

Figures

Equations

$${\hat{H}}_S=\sum _{i,\alpha }{u}_{\alpha }^i{\hat{S}}_{\alpha }^i+\sum _{i\ne j,\alpha ,\beta }{w}_{\alpha ,\beta }^{i,j}{\hat{S}}_{\alpha }^i{\hat{S}}_{\beta }^j,$$

$$\frac{d}{dt}⟨{\hat{S}}_{\alpha }^i⟩=\mathrm{i}⟨[{\hat{H}}_S,{\hat{S}}_{\alpha }^i]⟩.$$

$$⟨\hat{O}⟩=\mathrm{Tr}[\hat{\rho }\hat{O}]=\int dxdpW(x,p)O_W(x,p).$$

$$⟨\hat{O}⟩(t)\approx \int d{x}_{0}d{p}_{0}W(x_0,p_0)O_W(x_{\mathrm{cl}}(t),p_{\mathrm{cl}}(t)).$$

$$W(S_x^i,S_y^i,S_z^i)=\frac{1}{\pi S}{\mathrm{e}}^{-[(S_x^i{)}^2+(S_y^i{)}^2]/S}\delta (S_z^i-S),$$

$$⟨\hat{O}⟩(t)\approx \int \prod _{i,\alpha }d{S}_{\alpha ,0}^{i}W(S_{\alpha ,0}^i)O_W(\{S_{\beta ,\mathrm{cl}}^j(t)\}),$$

$$⟨{\hat{S}}_x^m{\hat{S}}_x^n⟩(t)$$

$$\ne ⟨{\hat{S}}_x^m⟩(t)⟨{\hat{S}}_x^n⟩(t),$$

$${\Lambda }_{\mu =1,\dots ,\mathcal{N}(\mathcal{N}-1)/2}^i$$

$$\text{for}\alpha >\beta ,1\le \alpha ,\beta \le \mathcal{N},$$

$${\Lambda }_{\mu =\mathcal{N}(\mathcal{N}-1)/2+1,\dots ,\mathcal{N}(\mathcal{N}-1)}^i$$

$$\text{for}\alpha >\beta ,1\le \alpha ,\beta \le \mathcal{N},$$

$${\Lambda }_{\mu =\mathcal{N}(\mathcal{N}-1)+1,\dots ,{\mathcal{N}}^2-1}^i$$

$$\sum _{\beta =1}^{\alpha }\left(|\beta ⟩{⟨\beta |}_i-\alpha |\alpha +1⟩{⟨\alpha +1|}_i\right)$$

$$\text{for}1\le \alpha <\mathcal{N},$$

$${\Lambda }_{\mathcal{D}}^i$$

$${\hat{O}}^i$$

$${\hat{H}}_{\Lambda }=\sum _{i,\mu }{u}_{\mu }^i{\Lambda }_{\mu }^i+\sum _{i,j\ne i,\mu ,\nu }{w}_{\mu ,\nu }^{i,j}{\Lambda }_{\mu }^i{\Lambda }_{\nu }^j,$$

$$\mathrm{i}\frac{d}{dt}⟨{\Lambda }_{\mu }^i⟩$$

$$=\sum _{\mu }{v}_{\nu }^i⟨[{\Lambda }_{\mu }^i,{\Lambda }_{\nu }^i]⟩+\sum _{\sigma ,j\ne i,\nu }{w}_{\sigma ,\nu }^{ij}⟨[{\Lambda }_{\mu }^i,{\Lambda }_{\sigma }^i{\Lambda }_{\nu }^j]⟩.$$

$${\lambda }_{\mu }^i\in \{a_{\mu }^i\}\text{with probabilities}$$

$$p_{\mu }^i({\lambda }_{\mu }^i=a_{\mu }^i)=\mathrm{Tr}\left({\hat{\rho }}^i|a_{\mu }^i⟩⟨a_{\mu }^i|\right).$$

$$⟨{\hat{O}}^m⟩$$

$$=\sum _{\mu }{c}_{\mu }^m\sum _{a_{\mu }^m}{p}_{\mu }^m({\lambda }_{\mu }^m=a_{\mu }^m){\lambda }_{\mu }^m=\sum _{\mu }{c}_{\mu }^m⟨⟨{\lambda }_{\mu }^m⟩⟩.$$

$$\hat{\mathcal{O}}$$

$$⟨\hat{\mathcal{O}}⟩(t)\approx \sum _{{\mu }_1,{\mu }_2,\dots ,{\mu }_N}{C}_{{\mu }_1,{\mu }_2,\dots ,{\mu }_N}⟨⟨\prod _i{\lambda }_{{\mu }_i}^i(t)⟩⟩.$$

$$⟨{\hat{O}}^m⟩(t)$$

$${\Lambda }_1$$

$${\Lambda }_3$$

$${\Lambda }_5$$

$${\Lambda }_7$$

$${\hat{H}}_{dd}$$

$${\xi }_R^2=\frac{2NS\min [\Delta S_{\perp }^2]}{|\vec{S}{|}^2},$$

$$\gamma (S,m_0)$$

$$f_{\pm }$$

$$\beta \approx -\frac{24E_0}{S^2(S+1{)}^2{\sum }_{i,j\ne i}{V}_{i,j}^2}.$$

$$n_m(\infty )=\frac{\mathrm{Tr}[{\hat{\rho }}_{\infty }{\hat{n}}_m]}{\Xi }.$$

$${\mathcal{C}}_S(t)$$

Full text

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 $S=1/2$ 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 $S>1/2$ 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 $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 Zaletel et al. (2015); Haegeman et al. (2016); Schachenmayer et al. (2013); Hauke and Tagliacozzo (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 $S$. The method is based on the well known truncated Wigner approximation (TWA) 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 $3$ spin-variables, but all $(2S+1)^{2}$ 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 $S=1/2$ our generalized discrete TWA approach (GDTWA), reduces to the previously proposed discrete TWA method (DTWA) Schachenmayer et al. (2015a), which has been remarkably succesful in predicting $S=1/2$ 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 $N$ spin-$S$ particles coupling to each other via two-body interactions can be written as

[TABLE]

with $\alpha,\beta=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, $\langle\hat{S}^{i}_{\alpha=x,y,z}\rangle$, can be obtained via the Heisenberg equations of motion (we set $\hbar\equiv 1$ throughout this paper):

[TABLE]

These equations generally depend on high order inter-spin correlations such as $\langle\hat{S}_{\alpha}^{i}\hat{S}^{j}_{\beta}\rangle$, $\langle\hat{S}_{\alpha}^{i}\hat{S}^{j}_{\beta}\hat{S}^{k}_{\gamma}\rangle$, $\dots$, for $i\neq j\neq 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, $\langle\hat{S}_{\alpha}^{i}\rangle$, e.g., $\langle\hat{S}_{\alpha}^{i}\hat{S}^{j}_{\beta}\rangle\approx\langle\hat{S}_{\alpha}^{i}\rangle\langle\hat{S}_{\beta}^{j}\rangle$. Then Eqs. (2) turn into $N$ coupled closed equations which can be easily solved numerically, and which correspond to equations for classical spin-variables111Note that more general models than Eq. (1) for $S>1/2$ can also depend on intra-spin terms such as quadratic fields $\propto(\hat{S}_{z}^{i})^{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.. 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 $\sim N^{2}$ and furthermore higher order corrections can typically lead to numerical instabilities. An alternative approach to simulate quantum correlations while retaining a complexity $\sim N$ 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 $\hat{O}$ 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 $\hat{O}$ can then be computed as an integral in phase space:

[TABLE]

Here, $W(x,p)$ is the Wigner function corresponding to the Weyl symbol of the density matrix $\hat{\rho}$: $W(x,p)=1/(2\pi)\int ds\,\bra{x+s/2}\hat{\rho}\ket{x-s/2}\exp(-{\rm i}ps)$, and represents a quasiprobability distribution for points in phase space. In the truncated Wigner approximation (TWA), the time evolution of $\langle\hat{O}\rangle$ is approximated from the initial Wigner function and the classical evolution of the phase space variables:

[TABLE]

Here, $x_{\rm cl}(t)$ and $p_{\rm cl}(t)$ are classical trajectories of $x$ and $p$, respectively, with $x_{\rm cl}(0)=x_{0}$ and $p_{\rm 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 $\hbar$ 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 $O_{W}(x_{\rm cl}(t),p_{\rm cl}(t))=x_{\rm cl}^{2}(t)p_{\rm cl}(t)$ corresponds to the observable $\hat{O}=(\hat{x}^{2}\hat{p}+2\hat{x}\hat{p}\hat{x}+\hat{p}\hat{x}^{2})/2$ Polkovnikov (2010).

For a system of $N$ coupled spin-$S$ particles, a natural way to formulate the TWA dynamics consists of replacing $\{\hat{x},\hat{p}\}\rightarrow\{\hat{S}^{i}_{\alpha}\}$, and using the $3N$ spin-variables $\{S^{i}_{\alpha}\}$ (with $\alpha=x,y,z$ and $i=1,2,\dots N$), corresponding to the Weyl symbols of $\{\hat{S}^{i}_{\alpha}\}$, as phase space variables. The trajectories $S^{i}_{\alpha,{\rm cl}}(t)$ correspond to those obtained from the mean-field approximation of Eq. (2) with $S_{\alpha,{\rm cl}}^{i}\equiv\langle\hat{S}^{i}_{\alpha}\rangle$. For particular initial states, the Wigner function can be easily found. Taking for example a product state $\prod_{i}\hat{\rho}^{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 form

[TABLE]

where $\delta(\cdot)$ 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 $\Delta(S_{x}^{i})^{2}=\Delta(S_{y}^{i})^{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 $\hat{S}^{i}_{x}$ and $\hat{S}^{i}_{y}$. For large $S\gg 1$, $\Delta S^{i}_{x}/S\to 0$, suggesting vanishing effects from quantum fluctuations.

For such initial states, the TWA evolution now reads

[TABLE]

where the Weyl symbol corresponds again to symmetrized observables. For example, for an inter-spin correlation such as $\hat{O}=\hat{S}_{x}^{n}\hat{S}_{x}^{m}$ at distinct sites $n\neq m$, simply $O_{W}(\{S_{\beta,\rm{cl}}^{j}(t)\})=S_{x,{\rm cl}}^{n}(t)S_{x,{\rm cl}}^{m}(t)$. 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 $m$ can now depend on the initial conditions of all other spin-variables via some function, e.g. $S_{x,{\rm cl}}^{m}(t)\equiv\mathcal{F}_{x}^{m}(\{S_{\beta,0}^{j}\})$. Therefore, for observables computed in the TWA, such as $\hat{S}^{m}_{x}\hat{S}^{n}_{x}$ with $m\neq n$,

[TABLE]

and correlations can build up with time. Generally, the factorization in the last equation would only hold if $\mathcal{F}_{x}^{m}$ would only depend on the variables of the spin $m$ itself.

II.2 Generalized discrete Truncated Wigner approximation (GDTWA)

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 $\propto(\hat{S}_{z}^{i})^{2}$ 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 $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 Bertlmann and Krammer (2008). This can be accomplished by noticing that for a spin-$S$ atom with $\mathcal{N}=2S+1$ spin states, its density matrix $\hat{\rho}_{i}$ consists of $\mathcal{D}=\mathcal{N}\times\mathcal{N}$ elements. Correspondingly, we can define $\mathcal{D}$ hermitian operators, $\Lambda_{\mu}$ (hats are dropped for clarity of notation), using the generalized Gell-Mann matrices (GGM) and the identity matrix $\mathbb{1}$222Note 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.:

[TABLE]

For each spin $i$, these matrices are orthonormal, ${\rm Tr}[\Lambda^{i}_{\mu}\Lambda^{i}_{\nu}]=\delta_{\mu,\nu}$, and constitute a complete local basis. Hence any single-spin operator can be represented via

[TABLE]

and $\mu=1,2,...,\mathcal{D}$. Consequently, more general Hamiltonians than Eq. (1), with one- and two-body terms can now be represented as

[TABLE]

with $\mu,\nu=1,2,\dots,\mathcal{D}$ and $i,j=1,2,\dots,N$. The evolution of expectation values of GGMs can now be computed as

[TABLE]

For each GGM $\Lambda^{i}_{\mu}$ we define a corresponding real phase space variable, $\lambda^{i}_{\mu}$. 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 motion follow from assuming a factorization between different sites $\langle\Lambda^{i}_{\mu}\Lambda^{j}_{\nu}\dots\Lambda^{k}_{\sigma}\rangle=\langle\Lambda^{i}_{\mu}\rangle\langle\Lambda^{j}_{\nu}\rangle\dots\langle\Lambda^{k}_{\sigma}\rangle$ for any nonequal site-index $i,j,\dots,k$ in Eq. (14), and by replacing $\lambda^{i}_{\mu}(t)\equiv\langle\Lambda_{\mu}^{i}\rangle$.

To define a phase space probability distribution for the inital state, we decompose each $\Lambda^{i}_{\mu}$ via its eigenvectors $\ket{a^{i}_{\mu}}$ with corresponding eigenvalues $a^{i}_{\mu}$, $\Lambda^{i}_{\mu}=\sum_{a^{i}_{\mu}}a^{i}_{\mu}\ket{a^{i}_{\mu}}\bra{a^{i}_{\mu}}$. As in spin-1/2 systems, where the Pauli matrices $\sigma_{x,y,z}$ can be measured to be $\pm 1$ in a projective measurement, the $a^{i}_{\mu}$ correspond to possible measurement outcomes of $\Lambda^{i}_{\mu}$. We will focus on the experimentally relevant case of initial product states $\hat{\rho}=\prod_{i}\hat{\rho}^{i}$. Then, for each $\hat{\rho}^{i}$ we can define a corresponding probability distribution for our phase space variables, which we limit to a discrete set of values

[TABLE]

The overall distribution factorizes for different variables on the same site and between sites, such that the probability for a configuration of all $\lambda^{i}_{\mu}$ being a certain combination of the eigenvalues $a_{\mu}^{i}$ for $i=1,2,\dots,N$ and $\mu=1,2,...,\mathcal{D}-1$ is given by $p(\{\lambda^{i}_{\mu}=a^{i}_{\mu}\})=\prod_{i,\mu}p_{\mu}^{i}(\lambda^{i}_{\mu}=a^{i}_{\mu})$.

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

[TABLE]

Here we defined the notation $\langle\!\langle\cdot\rangle\!\rangle\equiv\sum_{\{a_{\mu}^{i}\}}p(\{\lambda^{i}_{\mu}=a^{i}_{\mu}\})(\cdot)$, denoting the statistical average over any combination of eigenvalues for the $\lambda_{\mu}^{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 spin-operators on a site $m$, $\hat{S}^{m}_{x,y,z}$ and $\mathbb{1}$ do not form a complete local basis for any operator. Thus, for example the observable $(\hat{S}_{z}^{m})^{4}$ cannot be expanded as linear superposition of the matrices $\{\hat{S}^{m}_{x,y,z},\mathbb{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 $(\hat{S}_{z}^{m})^{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 $\lambda^{i}_{\mu}(0)\in\{a_{\mu}^{i}\}$ for all $\mu$ 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\bm{\lambda}(t)/dt=\mathcal{M}(\bm{\lambda}(t))\bm{\lambda}(t)$. Here $\bm{\lambda}$ is a vector of all $N\times(\mathcal{D}-1)$ phase-space variables $\lambda_{\mu}^{i}$, and the matrix $\mathcal{M}(\bm{\lambda}(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, $\lambda^{m}_{\nu}(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 $\Lambda_{\nu}^{m}$ at time $t$ as $\langle\Lambda_{\nu}^{m}\rangle(t)\approx\langle\!\langle\lambda^{m}_{\nu}(t)\rangle\!\rangle$.

If the problem is linear, i.e. if $\mathcal{M}(\bm{\lambda}(t))=\mathcal{M}$ 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 $\langle\Lambda_{\nu}^{m}\rangle(t)=\exp(\mathcal{M}t)\langle\!\langle\lambda^{m}_{\nu}(0)\rangle\!\rangle$, and since $\langle\!\langle\lambda^{m}_{\nu}(0)\rangle\!\rangle=\langle\Lambda^{m}_{\nu}\rangle(0)$ 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 $\langle\!\langle\lambda^{m}_{\nu}(t)\rangle\!\rangle$ will also depend on intra-spin correlations in the initial state, also within a single spin, e.g. $\langle\!\langle\lambda^{i}_{\mu}(0)\lambda_{\nu}^{i}(0)\rangle\!\rangle$, 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. $\langle\!\langle{\lambda^{i}_{\mu}(0)\lambda_{\nu}^{i}}(0)\rangle\!\rangle=\langle\!\langle{\lambda^{i}_{\mu}(0)}\rangle\!\rangle\langle\!\langle{\lambda_{\nu}^{i}(0)}\rangle\!\rangle$ with $\mu\neq\nu$. This is not necessarily true for the quantum mechanical intra-spin correlations of a generic single-spin state. For example, for some states $\hat{\rho}^{i}_{g}$, ${\rm Tr}[\hat{\rho}^{i}_{g}(\Lambda^{i}_{\mu}\Lambda^{i}_{\nu}+\Lambda^{i}_{\nu}\Lambda^{i}_{\mu})/2]\neq{\rm Tr}(\hat{\rho}^{i}_{g}\Lambda^{i}_{\mu}){\rm Tr}(\hat{\rho}^{i}_{g}\Lambda^{i}_{\mu})$, 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 $\hat{\rho}=\prod_{i}\ket{\alpha_{0}}_{i}\bra{\alpha_{0}}$ where $\alpha_{0}=1,2,\dots,\mathcal{N}$, and $\ket{\alpha_{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 $\hat{\rho}=\prod_{i}\ket{\alpha_{0}}_{i}\bra{\alpha_{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 $\lambda_{\mu}^{i}$ Wurtz et al. (2018); Davidson and Polkovnikov (2015), i.e. for all $\mathcal{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 Schachenmayer et al. (2015a).

To compute the time-dependent expectation value of an arbitrary multi-spin observable, $\hat{\mathcal{O}}$, we expand

[TABLE]

and approximate

[TABLE]

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, $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, $\hat{O}^{m}$. After selecting a single random configuration for all phase space variables $\{\lambda_{\mu}^{i}\}$ according to $p_{\mu}^{i}(\{\lambda_{\mu}^{i}=a_{\mu}^{i}\})$, we denote the subsequent classical evolution of the variables at site $m$ for this sample “$s$” as $\lambda_{\mu,{\rm cl}}^{m,[s]}(t)$. Then, the evolution of $\hat{O}^{m}$ is computed as

[TABLE]

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 in Eq. (19) 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 $\Lambda_{\mu}$ for each spin are proportional to standard Pauli matrices, $\Lambda_{x,y,z}=\sigma_{x,y,z}/\sqrt{2}$, each of which has eigenvalues $\pm 1/\sqrt{2}$, and eigenstates described by states aligned (anti-aligned) along the corresponding direction: $\ket{\uparrow_{x,y,z}(\downarrow_{x,y,z})}$. For a spin initially polarized as $\ket{\uparrow_{z}}$, the initial values of $\lambda_{x}$ in phase space take its two eigenvalues, $\pm 1/\sqrt{2}$, with equal probability ${\rm Tr}[\ket{\uparrow_{z}}\braket{\uparrow_{z}}{\uparrow_{x}}\bra{\uparrow_{x}}]=|\braket{\uparrow_{z}}{\uparrow_{x}}|^{2}=|\braket{\uparrow_{z}}{\downarrow_{x}}|^{2}=0.5$. In the same manner probabilities for $\lambda_{y}=\pm 1/\sqrt{2}$ are $0.5$, while $\lambda_{z}=1/\sqrt{2}$ is fixed (probability 1). Our GDTWA prescription thus reduces to the DTWA sampling introduced in Ref. Schachenmayer et al. (2015a).

For $S=1$, there are three states $\ket{+1},\ket{0},\ket{-1}$ for each spin, and $\Lambda_{\mu}$ consists of the identity matrix and 8 nontrivial Gell-Mann matrices:

[TABLE]

An initially polarized state $\ket{+1}$, can now be represented via probabilities for eigenvalues of $\Lambda_{1,2,...8}$. Specifically, with $1/2$ probability one chooses $\lambda_{1},\lambda_{2},\lambda_{4},\lambda_{5}$ at value $1/\sqrt{2}$, and with $1/2$ probability one chooses the value $-1/\sqrt{2}$. In the following, we denote this as $\lambda_{1,2,4,5}=\{-1/\sqrt{2},1/\sqrt{2}\}$, with probability $p_{1,2,4,5}=\{0.5,0.5\}$. Then, further we have fixed $\lambda_{3,6}=0$ with $p_{3,6}=1$, $\lambda_{7}=1/\sqrt{2}$ with $p_{7}=1$, and $\lambda_{8}=1/\sqrt{6}$ with $p_{8}=1$ (See Fig. 1 for an illustration of this distribution).

While the initial states discussed above are still spin coherent states, the state $\ket{0}$, 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 $\lambda_{1,3,4,6}=\{-1/\sqrt{2},1/\sqrt{2}\}$, $p_{1,3,4,6}=\{0.5,0.5\}$, $\lambda_{2,5}=0$, $p_{2,5}=1$, $\lambda_{7}=-1/\sqrt{2}$, $p_{7}=1$, and $\lambda_{8}=1/\sqrt{6}$, $p_{8}=1$.

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:

[TABLE]

where $V_{i,j}=V_{dd}(1-3\cos^{2}\theta_{i,j})/r_{i,j}^{3}$ is the strength of the interaction between two atoms $i$ and $j$, whose distance is $r_{i,j}$ and $\theta_{i,j}$ is the angle between the vector ${\bf 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}=\mu_{0}(\mu_{B}g)^{2}/4\pi$, with $\mu_{0}$ the magnetic permeability of vacuum, $\mu_{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 Yan et al. (2013); Lepoutre et al. (2019) $\ket{\psi(t=0)}=\otimes_{j}\exp({{\rm i}\theta\hat{S}_{y}^{j}})\ket{-S}_{j}$, with $\theta\in[0,\pi]$ the tipping angle, or ii) all atoms in the same Zeeman level de Paz et al. (2016a, 2013, b); Patscheider et al. (2019) $\ket{\psi(t=0)}=\otimes_{j}\ket{m_{0}}_{j}$, where $m_{0}=-S,-S+1,\dots,S,$ and $j=1,2,\dots,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)\sum_{j}|\braket{\psi(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 using 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 $\langle\hat{S}^{i}_{x,y,z}\rangle$, 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 $\langle\hat{S}_{x,y}^{n}\rangle$ 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, $\mathcal{S}^{c}_{\alpha}=\frac{1}{1-\alpha}{\log_{2}}{\rm Tr}[\hat{\rho}_{c}^{\alpha}]$, where $\hat{\rho}_{c}$ is the reduced density matrix of a single spin $c$, and $\alpha\geq 0$. In Fig. 3, we compare the evolution of the second Renyi entropy ($\alpha=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 $\theta$ 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 ($\theta\to 0$). This is an interesting observation that needs further investigation given that the $\theta\to 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 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),

[TABLE]

where $\Delta S_{\perp}^{2}=\langle\hat{S}_{\perp}^{2}\rangle-\langle\hat{S}_{\perp}\rangle^{2}$ is the variance measured perpendicular to the collective spin $\vec{S}$ direction, and $|\vec{S}|=\sqrt{\langle\hat{S}^{x}\rangle^{2}+\langle\hat{S}^{y}\rangle^{2}+\langle\hat{S}^{z}\rangle^{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, $\Delta S_{\perp}^{2}=NS/2$, $\xi_{R}^{2}=1$. On the other hand $\xi_{R}^{2}<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 spin-observables such as $\langle(\hat{S}^{x})^{2}\rangle$. In various benchmark situations we find that this modified procedure provides much more precise predictions for the spin-squeezing parameter 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.

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 $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. Patscheider et al. (2019), the effective dipolar exchange strength given by $V\gamma(S,m_{0})$ can be used as the characteristic energy scale of the spin dynamics, with $\gamma(S,m_{0})$ given by

[TABLE]

where $m_{0}$ is the initially populated spin level. Correspondingly we define $\tilde{V}=V\gamma(S,m_{0})$. As shown in Fig. 5(b), for various values of $S$, the initial dynamics happens at a characteristic rate $1/\tilde{V}$.

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 $\sim 200$ atoms, which corresponds to a full Hilbert space $\sim 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 non-integrable system far away from a localized phase 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, $\hat{\rho}_{\infty}={\rm exp}(-\beta\hat{H}_{T}-\mu\hat{S}^{z})/\Xi$, with $\Xi={\rm Tr}[{\rm exp}(-\beta\hat{H}_{T}-\mu\hat{S}^{z})]$, $\hat{H}_{T}=\hat{H}_{dd}$, and $\hat{S}^{z}=\sum_{i}\hat{S}_{z}^{i}$. The effective inverse temperature $\beta$ and chemical potential $\mu$ 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}=\sum_{i,j\neq i}V_{i,j}m_{0}^{2}/2$. In a high temperature limit $\beta\to 0$, energy conservation leads to

[TABLE]

For a 3D lattice where dipolar interactions vary in sign, $E_{0}\approx 0$ and thus we take $\beta=0$ in the following. Then $\mu$ can be easily determined from $\sum_{m=-S}^{S}m{\rm exp}(-\mu m)/\Xi=m_{0}$, with $\Xi=\sum_{m=-S}^{S}{\rm exp}(-\mu m)$, and we can obtain the equilibrium population in a spin $m$ state:

[TABLE]

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 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 $S=1/2$ Acevedo et al. (2017) and $S=3$ Lepoutre et al. (2019) models. Note that while our simulations do not constitute a general proof of thermalization in the GDTWA, our results here suggest 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 Távora et al. (2016, 2017); 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 $n_{m_{0}}(t)-n_{m_{0}}(\infty)\sim 1/t^{2}$, with $n_{m_{0}}(t)$ the population on state $m_{0}$ at time $t$.

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 $\mathcal{S}^{c}_{2}=-\log_{2}{\rm Tr}[\hat{\rho}_{c}^{2}]$ for a subsystem of a single spin. The entropy quantifies 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 $\hat{\rho}_{c}$ chosen as the reduced density matrix of a central spin in the cube. Since our initial state is a pure state with $\mathcal{S}_{2}^{c}=0$, an increasing $\mathcal{S}_{2}^{c}(t)>0$ indicates a build-up of entanglement in the system. For all values of $S$, $\mathcal{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 can feature richer entanglement-related quantum effects. Their increased local Hilbert space size allows larger entanglement entropies, and in our setup we find a quick approach to the maximally allowed Renyi entropy. Since for dipolar systems, such as the magnetic atoms discussed above, larger $S$ is associated with a larger net interaction strength, a larger spin also implies a faster generation of entanglement in the spin dynamics as confirmed in Fig. 6.

IV.3 Dynamics of contrast and intra-spin correlations

Another useful probe of many-body effects is the contrast of spin coherence, defined as $\mathcal{C}_{S}=\sqrt{\langle\hat{S}^{x}\rangle^{2}+\langle\hat{S}^{y}\rangle^{2}}/(NS)$, with $\hat{S}^{x,y,z}=\sum_{j}\hat{S}_{j}^{x,y,z}$. 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 $\mathcal{C}_{S}$ under the dipolar interactions in Hamiltonian (21), starting from the type i) spin coherent initial state with tipping angle $\theta=\pi/2$. Using a short-time expansion analysis, we find that the contrast initially decays as Hazzard et al. (2014b)

[TABLE]

Compared with the dynamics discussed in the previous section, where a larger spin $S$ can provide an enhanced dipole-dipole interaction $\propto\tilde{V}\propto S^{2}V$, here the contrast decays slower than the timescale $\propto 1/\tilde{V}$. Instead, the characteristic interaction strength for the contrast dynamics is $\tilde{V}_{c}=\sqrt{S}V$. This is seen in Fig. 7(a), where as a function of $\tilde{V}_{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 Vitagliano et al. (2014). As discussed in Sec. II, since in a Gaussian TWA only the spin operators $\hat{S}_{x,y,z}^{j}$ are involved, such intra-spin correlations generally cannot be well captured. In the initial state, the GDTWA reproduces the exact results, differing from those of the Gaussian sampling, as also discussed in B. Furthermore, we find that the discrepancy between the Gaussian TWA and the GDTWA increases with time, as demonstrated in Fig. 7(b). There, to illustrate this effect, we compare the evolution of $\langle\sum_{j}(\hat{S}_{z}^{j})^{n}\rangle$ ($n=4,6$), computed with GDTWA and Gaussian TWA for three spin-variables. The Gaussian TWA provides incorrect initial values, and this error propagates over time.

IV.4 Quantum spin squeezing

Lastly, to further characterize the quantum correlations generated during the dipolar spin dynamics, we study the spin squeezing parameter $\xi_{R}^{2}$ introduced in Sec. III for ensembles of spins with various $S$. It is straightforward to show that $\xi_{R}^{2}\geq 1/(2NS)$. Hence, to make a fair comparison for systems of different spin $S$, in Fig. 8, we keep $NS$ constant and calculate $\xi_{R}^{2}$ 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 Vitagliano et al. (2014), 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({\mathcal{N}})$ 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.

A general drawback of the TWA is the lack of a clear convergence parameter, and consequently the difficulty to assess the validity at long times. While this validity can be verified by comparison in analytical limits, with exact calculations for feasible system sizes, or by corroborating it with experimental data, it is desirable to have a more systematical validation procedure. With our GDTWA approach such a procedure could be achieved by systematically grouping multiple spin-$S$ particles (local Hilbert size of $\mathcal{N}=2S+1$) into pairs of effectively larger spin-$\tilde{S}$ particles (local Hilbert space size $\mathcal{N}^{2}$) or larger clusters of $s_{c}$ spin-$S$ particles (local Hilbert space size $\mathcal{N}^{s_{c}}$). In such a way the exact quantum dynamics is preserved within the clusters, while only inter-cluster correlations are subjected to the TWA approximation. This comes at the cost of changing the number of classical equations from $N\times(\mathcal{N}^{2}-1)$ to $N_{c}\times(\mathcal{N}^{2s_{c}}-1)$, where $N_{c}$ is the number of clusters. By comparing the dynamics obtained for increasing cluster sizes $s_{c}$ (the simulation would hypothetically become exact in the extreme limit $N=s_{c}$) and different cluster choices, one could benchmark the convergence of the phase space method in regimes where no other methods are available for verification. For example, such an approach was recently suggested and analyzed 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 $\langle\!\langle\lambda^{i}_{\mu}(0)\lambda_{\nu}^{j}(0)\rangle\!\rangle$, $\langle\!\langle\lambda^{i}_{\mu}(0)\lambda_{\nu}^{j}(0)\lambda_{\kappa}^{k}(0)\rangle\!\rangle$ 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 $\ket{\alpha_{0}}$ with $\alpha_{0}=1,2,\dots,\mathcal{N}$, we mean any eigenstate of the diagonal GGMs defined in Eqs. (10), which are eigenstates of $\hat{S}_{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 $\hat{\rho}=\ket{\alpha_{0}}\bra{\alpha_{0}}$,

[TABLE]

where we used the expansion of the GGM into the orthonormal eigenbasis $\Lambda_{\mu}=\sum_{a_{\mu}}a_{\mu}\ket{a_{\mu}}\bra{a_{\mu}}$. Since the discrete distributions defined in Eq. (15) are furthermore independent for different $\mu$, correlations between different GGMs always factorize:

[TABLE]

with arbitrary powers $n_{1},n_{2},\dots,n_{\mathcal{D}}\geq 0$. We will show that this property also holds for correlations for the diagonal quantum state.

We proceed by showing that for $\mu\neq\nu$

[TABLE]

The symmetrization super-operator, $\mathcal{S}$, is necessary since generally different GGMs cannot be measured simultaneously as they do not always commute. Note that (29) is always true if $\ket{\alpha_{0}}$ is an eigenstate of $\Lambda_{\mu}$ and $\Lambda_{\nu}$. 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 $\ket{\alpha_{0}}$ there are only $2(\mathcal{N}-1)$ GGMs for which $\ket{\alpha_{0}}$ is not an eigenstate,

[TABLE]

with $\xi=1,2,\dots,\alpha_{0}-1,\alpha_{0}+1,\dots,\mathcal{N}$. For all other off-diagonal GGMs, $\ket{\alpha_{0}}$ is an eigenstate with zero eigenvalue. We therefore have to only show that Eq. (29) also holds if at least one of the matrices $\Lambda_{\mu}$ and $\Lambda_{\nu}$ is from Eqs. (A). It is straightforward to compute that for $\mu\neq\nu$

[TABLE]

where $\tilde{\Lambda}_{\mu}$ is any GGM not from (A). Therefore, (29) 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 $\ket{\alpha_{0}}$ as eigenstate, or it is proportional to itself, therefore also ($\mu\neq\nu$)

[TABLE]

holds for arbitrary integer powers of $n_{\mu}$ and $n_{\nu}$.

Lastly, it’s straightforward to see that the pairwise property (35) also leads to a factorization of any other symmetrized product of multiple GGMs. To realize this, consider that

[TABLE]

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, $\Lambda_{\xi}$. If $\Lambda_{\xi}$ commutes with $(\Lambda_{\kappa})^{n_{\kappa}}$ we do not need to symmetrize, otherwise $\kappa\neq\xi$ and we can again use Eq. (35) to see that $\bra{\alpha_{0}}\mathcal{S}\Lambda_{\xi}(\Lambda_{\kappa})^{n_{\kappa}}\ket{\alpha_{0}}$ factorizes. From a repeated application of this argument it follows that

[TABLE]

In particular, for states $\ket{\psi}$ that are not diagonal, but can be obtained from a diagonal state via any single-site unitary transformation ($\hat{U}^{\dagger}\hat{U}=\mathbb{1}$), $\ket{\psi}=\hat{U}\ket{\alpha}$, we can transform the Hamiltonian into the rotated basis $\hat{H}^{\prime}=\hat{U}^{\dagger}\hat{H}\hat{U}$, and then obtain the equations of motion for GGMs in a similar way as Eq. (14), with $\hat{H}$ replaced by $\hat{H}^{\prime}$. The dynamics of an observable of interest $\hat{O}$ can the be found from evaluating $\hat{O}^{\prime}=\hat{U}^{\dagger}\hat{O}\hat{U}$. Since in this rotated basis the initial state becomes diagonal, the expectation values of $\hat{O}^{\prime}$ can be correctly described with the discrete sampling, and thus $\bra{\psi}\hat{O}\ket{\psi}$ in the original basis can be obtained. For the type i) initial states (used in the main text) with arbitrary $\theta$, we can use $\hat{U}=e^{i\theta\hat{S}^{y}}$ and perform sampling in the rotated basis with initial state $\ket{-S}$. For initial states that are a coherent superposition of several Zeeman levels, $\ket{\psi}=\sum_{m}c_{m}\ket{m}$, with $m=-S,-S+1,...S$, $c_{m}$ real numbers and $\sum_{m}c_{m}^{2}=1$, the unitary transformation can be constructed from $\hat{U}=\prod_{m=-S+1}^{S}{\rm e}^{-{\rm i}{\rm arcsin}[c_{m}/\sqrt{1-\sum_{m^{\prime}=S}^{m+1}c_{m^{\prime}}^{2}}]\hat{A}_{-S,m}}$, with $\hat{A}_{-S,m}=\ket{-S}\bra{m}+\ket{m}\bra{-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: $\hat{U}_{N}=\otimes_{j}\hat{U}_{j}$ with $\hat{U}^{\dagger}_{j}\hat{U}_{j}=\mathbb{1}_{j}$.

Appendix B Failure of state-tomography for Gaussian states

Let’s take a Gaussian state polarized along the $x$-axis, following from a $\pi/2$ rotation of the state from Eq. (5) (site index dropped):

[TABLE]

This Gaussian distribution only reproduces the mean and the variance of the spin-operators correctly, but not the higher moments $\overline{S_{z}^{n}}=\int_{-\infty}^{\infty}d\bm{S}\,W(S_{x},S_{y},S_{z})S_{z}^{n},$ with $n>3$. For example, the first $4$ moments of the $z$ component are $\overline{S_{z}^{0}}=1,\overline{S_{z}^{1}}=0,\overline{S_{z}^{2}}={S}/{2},\overline{S_{z}^{3}}=0,\overline{S_{z}^{4}}={3}S^{2}/4$. In contrast, the exact quantum mechanical spin coherent state is

[TABLE]

and gives rise to moments distributed according to a binomial distribution. In particular, the exact $4$-th moment is $\bra{\rightarrow}\hat{S}_{z}^{4}\ket{\rightarrow}={3}S^{2}/4-{S}/{4}$ and differs from the Gaussian one by $S/4$.

The differences in the higher order moments has profound consequences when performing a state-tomography for large spin density matrices. For example, the diagonal density matrix elements (i.e. the populations) can be expressed by powers of the spin-$z$ matrix. In the following let’s take as an example the “central” Zeeman state $m_{c}=0$ ($S$ even) or $m_{c}=+1/2$ ($S$ odd). Then it is straightforward to see that

[TABLE]

We see that for $S\geq 2$ the 4th order moment starts to play a role if we want to evaluate the expectation values for the $\hat{S}_{z}^{n}$ operators from a probability distribution. In particular, assuming that we choose a Gaussian distribution, then for $S=2$

[TABLE]

while the exact population of the state $\ket{m_{c}}$ is

[TABLE]

This difference is significant and it does not vanish with increasing spin. Fig. 9 demonstrates the failure of the Gaussian sampling for obtaining Zeeman level populations for various $S\geq 2$. Especially, the large integer spin cases are very badly re-produced by a Gaussian distribution. In general, this shows that already the diagonal elements are not faithfully reproduced from a Gaussian sampling for an initial state, and the same arguments also hold for off-diagonal density-matrix terms, which in the case of $S\geq 2$ will require expectation values of $\hat{S}_{x}^{\geq 4},\hat{S}_{y}^{\geq 4}$. Note that in contrast, our discrete distribution for diagonal states (including spin coherent states) is exact and provides correct predictions for moments of all orders (see A).

Consequently, also an attempt to simulate the time-evolution of Zeeman level populations using the state-tomography from expectation values of powers of $\hat{S}_{j}^{z}$ gives wrong results for $S\geq 2$ as demonstrated in Fig. 10. There we perform the same simulation as in Fig. 2(a). The Gaussian sampling for the three spin-variables can reproduce results for dynamics of the spin-components very well (shown here: $\sum_{j}\langle\hat{S}_{x}^{j}\rangle/N$). In contrast, here the time-dependence of the populations has to be constructed from the time-dependent evolution of observables corresponding to the powers $(\hat{S}_{j}^{z})^{\eta}$ via expansions such as (43). While already the initial value for the powers with $\eta>3$ are badly reproduced (Fig. 9), the sub-sequent time-evolution then gives rise to even unphysical evolution of populations, reaching negative values or values larger than one already at short times. Note that the inital state of Fig. 2(b) cannot be described by a Gaussian distribution.

Appendix C Renyi entropy for smaller tipping angles ($\theta<\pi/2$)

In Fig. 3(a) we have shown that for $S=3$ and tipping angle $\theta=\pi/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 $\theta=0.2\pi$, these comparisons show a good agreement between the GDTWA and ED, suggesting that generally for $\theta$ not too small, the GDTWA approach works well for calculating the second Renyi entropy. Surprisingly, in the small $\theta$ 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 $\hbar$). Indeed this is corroborated by Fig. 11(b), where we show the case of $\theta=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 $\propto(Vt)^{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 $\xi_{R}^{2}$, we need to evaluate the expectation value of linear and quadratic operators such as $\hat{S}^{x}$ and $(\hat{S}^{x})^{2}$. In the GDTWA, we first evaluate the expansion coefficients for the general operators in the GGM basis [Eq. (17)]:

[TABLE]

Then, we approximate the quantum expectation values from averages over the classical trajectories: [Eq. (18)]:

[TABLE]

Alternatively, in the expansion of $(\hat{S}^{x})^{2}$ in Eq. (47) one could also ignore the distinct treatment given to the diagonal terms, and instead compute the expectation value $\langle(\hat{S}^{x})^{2}\rangle$ as

[TABLE]

In principle the first term in Eq. (49) takes the single-spin part of the collective observable into account more precisely. Note, e.g. that for $S=1/2$ it correctly reproduces the constant value for $\langle\hat{\sigma}_{x}^{2}\rangle=\langle\mathbb{1}\rangle=1$. Nevertheless, by comparisons with ED calculations, we find that in practice, especially for finite range interactions, Eq. (50) describes the evolution of $\xi_{R}^{2}$ quantitatively better than Eq. (49). This is illustrated in Fig. 12, where a comparison with ED shows a nearly perfect agreement using Eq. (50), while a clear deviation on the logarithmic scale is visible when using Eq. (49). Therefore, in the main text we adopted Eq. (50) for evaluating spin squeezing.