Dynamics of correlations in two-dimensional quantum spin models with long-range interactions: A phase-space Monte-Carlo study

Interacting quantum spin models are remarkably useful for describing different types of physical, chemical, and biological systems. Significant understanding of their equilibrium properties has been achieved to date, especially for the case of spin models with short-range couplings. However, progress towards the development of a comparable understanding in long-range interacting models, in particular out-of-equilibrium, remains limited. In a recent work, we proposed a semiclassical numerical method to study spin models, the discrete truncated Wigner approximation (DTWA), and demonstrated its capability to correctly capture the dynamics of one- and two-point correlations in one dimensional (1D) systems. Here we go one step forward and use the DTWA method to study the dynamics of correlations in 2D systems with many spins and different types of long-range couplings, in regimes where other numerical methods are generally unreliable. We compute spatial and time-dependent correlations for spin-couplings that decay with distance as a power-law and determine the velocity at which correlations propagate through the system. Sharp changes in the behavior of those velocities are found as a function of the power-law decay exponent. Our predictions are relevant for a broad range of systems including solid state materials, atom-photon systems and ultracold gases of polar molecules, trapped ions, Rydberg, and magnetic atoms. We validate the DTWA predictions for small 2D systems and 1D systems, but ultimately, in the spirt of quantum simulation, experiments will be needed to confirm our predictions for large 2D systems.


Introduction
An important advance towards understanding non-equilibrium phenomena has been made possible by recent advances in cooling, trapping and manipulating atomic, molecular and optical (AMO) systems [1]. In contrast to solid state systems, where studying equilibrium situations is the default approach (given their complex environment and fast relaxation rates), AMO systems provide a unique platform to observe and investigate non-equilibrium dynamics in strongly interacting many-body systems [2]. Their high dynamic tunability even during the course of an experiment, their decoupling from the external environment and their characteristic low-energy scales, lead to long non-equilibrium time-scales over which the system can be followed almost in real time. Quantum quenches, i.e. the dynamics induced by abruptly changing parameters of the system, is currently a common protocol used to probe AMO systems.
A complementary capability, which is emerging as one of the most promising opportunities offered by modern AMO physics is the ability to engineer interatomic interactions different from the standard contact and isotropic interactions arising from ultracold s-wave collisions. At the heart of this capability are recent experimental developments on controlling AMO systems with complex internal structure and with enlarged sets of degrees of freedom such as polar molecules [3], trapped ions [4][5][6], magnetic atoms [7][8][9], Rydberg atoms [10][11][12][13][14][15], and alkaline earth atoms [16][17][18][19][20]. All these systems have in common that they can exhibit long-range interactions. This experimental progress is opening new frontiers, and at the same time demanding for improved theoretical techniques that are capable of dealing with the complicated nonequilibrium dynamics of long-range interacting systems.
In a prior work [21], we proposed a semiclassical phase-space method to study nonequilibrium quantum dynamics. We used this numerical method, that we named the discrete Truncated Wigner Approximation (dTWA), to study the dynamics of single particle observables and correlation functions after a quench. We benchmarked the dTWA in one dimensional spin models with numerically exact time-dependent density matrix renormalization group calculations (t-DMRG) [22][23][24] and found excellent agreement.
In this work, we generalize the calculations to 2D Ising and XY spin models with various ranges of interactions. We study the time-evolution in a setup that is equivalent to a Ramsey-type procedure as realized in recent experiments. This dynamical protocol has been used, for example, to observe dipolar spin-exchange interactions in ultracold molecules [3,25,26], to benchmark the Ising dynamics of hundreds of trapped ions [5], and to precisely measure atomic transitions as well as many-body interactions in optical lattice clocks [17,18,27,28]. Using the dTWA we compute the dynamics of the collective spin as well as spatially resolved two-point correlation functions. Since the applicability of the t-DMRG method becomes limited in 2D, to benchmark the dTWA we perform numerical comparisons with small systems (where exact diagonalization is possible) and with the analytically solvable Ising case [29][30][31]. We then extend the calculations to large spin systems, and find remarkably sharp changes in behavior at a critical range of interactions.
The dynamics of the two-point correlations is directly linked to the speed of propagation of information in quantum many-body systems, a topic at the heart of quantum information science and currently subjected to intensive investigation. For systems with only short-range interactions, there is a well understood bound (derived by Lieb and Robinson) that limits correlations to remain within a linear effective "light cone" region [32,33]. On the contrary there are many open questions about what limits the propagation of information in quantum many-body systems with long-range interactions [34][35][36][37][38][39][40][41]. Here we use the dTWA method to determine the shape of the causal region and the speed at which correlations propagate. Our calculations and natural generalizations of them (e.g. local instead of global quenches) should be testable in experiments with polar molecules and trapped ions in the near future.
This paper is organized as follows: In section 2 we introduce the models that we study. In section 3 we review the dTWA technique that was introduced in Ref. [21]. In section 4 we benchmark the dTWA method by comparing the dynamics of singlespin observables and correlation functions with exact solutions. In section 5 those calculations are extended to large systems where currently no other method is applicable. Finally in section 6 we conclude and provide an outlook.

Spin models and dynamics
We will focus our attention on Hamiltonians that fall under the generic heading of spin-1/2 XXZ models given by ( = 1) where the sum extends over all pairs of sites of an arbitrary lattice,σ x,y,z i are Pauli matrices for the spin on site i, J ij = J ji and J ii = 0. In our analysis the interactions are assumed to decay as a function of the distance with a decay exponent α, J ⊥,z ij ≡ J(a/|r ij |) α . Here, r ij is the vector connecting spins on sites i and j and a is the lattice spacing. We concentrate our study on two specific cases, Ising (J ⊥ ij = 0) and XY (J z ij = 0) interactions. We consider a general 2D grid, i.e. a lattice with N x × N y = M sites with spins at positions r i = a(n x , n y ) with integer n x,y . The lattice spacing a is set to 1 throughout this paper.
Spin-1/2 XXZ models broadly describe a variety of physical systems. For instance, in the AMO context, XXZ spin Hamiltonians have been used to model the dynamics of ultracold molecules in optical lattices, trapped ions, Rydberg atoms, neutral atoms in optical clocks, and ultracold magnetic atoms. A summary of how these models are realized in those systems can be found in Ref. [42]. Here we describe the two most relevant ones for this work: ultracold polar molecules and trapped ions.
In ultracold polar molecules pinned in optical lattices, the spin-1/2 degree of freedom can be encoded in two rotational states, and the spin-spin couplings are generated by dipolar interactions. The difference in dipole moments between the two states (which arises in the presence of an electric field) generates the Ising term while transition dipole moments between the two rotational states (which can exist even in the absence of an electric field) give rise to the spin-exchange terms [43]. The ratio between the Ising and XY couplings can be manipulated using electromagnetic fields [44][45][46][47][48][49]. The case without electric field which implements the pure XY model has been recently realized with KRb polar molecules in a 3D optical lattice [3,26]. In general, dipolar interactions are long-ranged and spatially anisotropic. In a 2D geometry, however, they become isotropic if the electric field that sets the quantization axis is set perpendicular to the plane containing the molecules. This is the case considered throughout this paper.
Crystals of 2D self-assembled trapped ions can also be used to implement specific cases of equation (1), cf. [5]. By addressing the ions confined by a Penning trap with a spin-dependent optical potential, the vibrations of the crystal mediate a long-range Ising interaction that can be approximately described by a power-law with 0 ≤ α < 3 [4,5,[50][51][52][53]. To engineer an XY model, one needs to add a strong transverse field that projects out the off-resonant terms in the Ising interactions that change the magnetization along the field quantization direction [39,40].
The dynamic procedure considered here, is identical to a Ramsey spectroscopy setup. It has been implemented in various recent experiments as a diagnostic tool for interactions [26]. It consists of preparing an initial state with all spins aligned (at time t = 0) along a specific direction, here we consider it to be the x direction, i.e. |ψ( Then this initial state evolves under the Ising or XY Hamiltonian (1) for a time t, leading to the state |ψ(t) = exp(−itĤ) |ψ(t = 0) . Afterwards one measures expectation values of an observable with the time-evolved state, ψ(t)|Ô |ψ(t) . In this paper, we focus on two-point correlations and the collective spin along x as observables.

The dTWA method
Phase space methods, such as the Truncated Wigner Approximation (TWA), solve the quantum dynamics approximately by replacing the time-evolution by a semi-classical evolution via classical trajectories. The quantum uncertainty in the initial state is accounted for by an average over different initial conditions [54,55], determined by the Wigner function. Although the TWA approximation was initially developed to deal with systems with continuous degrees of freedom, it has also been adopted to treat collective spin models. In this case the standard method has been to approximate the Wigner function by a continuous Gaussian distribution that facilitates the sampling of trajectories [55]. This continuous approximation, however, misses important aspects inherent to the discrete nature of spin variables and it is inapplicable to systems with finite-range interactions. To deal with more generic types of spin models, recently we proposed instead to sample a discrete Wigner function for each spin and named this approach the dTWA method [21]. In this section we present an overview of the dTWA method. For details the reader is referred to Ref. [21].
Operators and wave functions on the Hilbert space of a quantum system can be, equivalently, represented (mapped) on a classical phase space. There, any operatorÔ corresponds to a real-valued function of the classical phase-space variables, O W (a socalled Weyl symbol). The phase-space function corresponding to the density matrix is precisely the Wigner function w. For continuous variables p, q in one dimension (for simplicity of presentation), the expectation value of an operator can be exactly represented as Ô (t) = dpdq w(p, q; t)O W (p, q). For quantum systems with discrete degrees of freedom, one can introduce a "discrete phase-space" in various ways, see [56] and references therein. Here we use the representation of Wootters [56,57]. For a single spin-1/2, it uses four distinct phase-points, and all phase-space functions are thus defined as 2×2 matrices. In our approximation, the phase-space of N spins factorizes into a product of N phase-spaces for each individual spin.
In both continuous and discrete variables however, it is not possible to compute the time-evolution exactly. The spirit of the TWA and dTWA is to take quantum fluctuations into account only to lowest order (in ) [55]. In particular, we switch to a "Heisenberg picture", such that the Wigner function does not evolve in time (i.e. it is fixed to the initial state) while the operator-functions are time-dependent. The approximation that we make is to assume that the operators in phase space follow their classical evolution: where γ runs over the points of the discrete phase space, w(γ; 0) is the Wigner function at t = 0 on the discrete many-body phase space, and O W,cl. (γ; t) is the classically evolved operator-function (Weyl symbol) that corresponds to our observable. Equation (2) is solved numerically by choosing a large number n t of random initial spin-configurations, with probability according to w(γ; 0). Each of this "Monte-Carlo trajectories" is evolved independently following the classical equations of motion (see below). The expectation value in equation (2) is calculated by averaging. We find that the number of required trajectories, n t , does not depend on the system size, but rather on the observable under consideration.
In practice, to apply the truncated Wigner approximation we have to compute the classical equations of motion for the spin components of each spin i: s x i , s y i , s z i . One way to do this is to replace spin operators σ by classical variables ( s = σ ) in the Hamiltonian and to compute the Poisson bracket [55], givingṡ δ i = 2 β δβγ s γ i ∂H ∂s β i with the fully antisymmetric tensor. Alternatively, the same classical (mean-field) equations of motion can be obtained from a product state ansatz of the density matrix [21]. For the Ising interaction Hamiltonian the classical equations for the spin components are given by:ṡ where we introduced the quantity β δ=x,y,z n ≡ m J z n,m s δ=x,y,z m which can be interpreted as an effective magnetic mean-field acting on spin n induced by the other spins. For the XY interaction the classical equations of motion are given by: In practice, applying the dTWA means to solve these equations of motion n t times, while each time choosing a different random initial configuration. For our particular initial state (pointing along the x direction), the prescription for the correct sampling is to randomly pick values of the orthogonal spin-components for each spin from s y n (0), s z n (0) ∈ {−1, +1}.

Contrast
Before discussing results obtained by the dTWA method for the propagation of correlations through a 2D system, we need to consider how well this approximate method performs in such a setting. As a starting point we first consider simpler single particle observables. One observable with immediate relevance to AMO experiments is the "contrast" or amplitude of the oscillations in a Ramsey experiment, which for our initial state is given byŜ x = iσ x i . In figure 1 we compare the time-evolution of Ŝ x to exact solutions. For Ising interactions, J ⊥ ij = 0 in equation (1), exact analytical expressions for the dynamics exists [29][30][31]. We can thus compare our dTWA results in a large 31 × 31 system. In contrast, for XY interactions [J z ij = 0 in equation (1)], no analytical solution is known in 2D. Therefore, in this case we have to resort to comparisons with a numerically exact diagonalization (ED) method, which is limited to small systems sizes. Here, we choose a 4 × 5 lattice.
To cover different regimes, in figure 1 we consider the Ising (panels a,c) and the XY (panels b,d) cases with two different power-law decay exponent, α = 3 (panels a, b), and α = 1 (panels c, d). Remarkably, the Ising dynamics is exactly covered by the dTWA approximation, an agreement that can be rigorously justified [21]. For the XY case we also find excellent agreement. While for α = 3 small numerical differences are visible, for α = 1, the different curves are nearly indistinguishable.

Spatial correlations
We now turn to checking the capability of the dTWA to describe the time evolution of spatial two-point correlations. We again first consider the case of Ising interactions where we can compare the dTWA to an exact analytical solution in large systems [29][30][31]. In particular we consider a 31 × 31 square lattice geometry. We study the time-dependence of connected correlation functions between sites n and m; we calculate the correlations from the central spin of the system, r n = (16, 16), along the y-direction, r m = (16, 16 + j) (note that results along the x-direction are identical). The comparisons for Ising interactions are summarized in figure 2. We find that in particular for substantially long-range interactions (case α = 1) the dTWA gives excellent results. Only for j ∼ 1, 2 there are small quantitative differences (the fact that for large j the error becomes small can be analytically understood [21]). For shorter range interactions (case α = 3) the discrepancies become larger. Nevertheless, the spreading of the correlations is still very well (quantitatively) reproduced by the dTWA. Note that the fluctuations around zero for the dTWA solution are statistical because of the finite number of trajectories, which is of order O(10 5 ) in all our calculations.
For the XY model we compare the dTWA against exact diagonalization results in a small 4 × 5 lattice. Due to the small system size, in order to observe any amount of linear spreading of correlations we have to calculate the correlations from the corner of the system. Specifically, in figure 3 we calculate the correlations between site n = i 0 with r i 0 = (1, 1) and sites j along the y-direction with coordinates r jy = (1, 1 + j). In this XY case we also find that the dTWA works impressively well, especially for j > 1.
Most importantly we find that although some oscillation seem not to be well reproduced, the dTWA accurately captures the spreading of the correlations and the shape of the "light-cone" boundary as seen in figure 3b. The "light-cone" boundary at a threshold value C thres is visualized in figure 3 by a contour plot. Physically it corresponds to the propagation time required to reach a correlation value of C thres between two spins separated by a distance j. Here we set C thres = 0.05. In contrast to the Ising case where the propagation of correlations barely follows a light-cone like spread with an expanding causal region, in the XY model the situation is more interesting. While for α = 1, the correlations build up throughout the system almost instantaneously, there is a clear change in behavior when going to shorter ranged interactions. In the case of α = 3 a clearly bounded causal region emerges [32,38], and is even visible in the small system calculation shown in figure 3b. However, in such a small system boundary effects are expected to play an important role. We point out that the speed of the propagation of correlations in figure 3 is essentially identical for C xx i,j and C yy i,j . However, the dTWA result for C xx i,j , shows slightly larger discrepancies from the exact solution than C yy i,j . We also checked the evolution of C zz i,j (in the XY case), which exhibits the same type of correlation spreading and is equally well reproduced by the dTWA. However, in our case this particular correlation is much smaller in magnitude than C yy i,j and features additional oscillations. Given the excellent predictions of the dTWA, especially for C yy i,j correlations, we now use the dTWA method to compute C yy i,j in larger systems, and study the maximum speed at which correlations can propagate in the XY model.

dTWA predictions for spreading of correlations in large systems
As in our previous Ising examples, we now focus on the large 31 × 31 square lattice geometry and study the evolution of C yy ic,jy . In figure 4 we summarize our large system results for decay exponents α = 1, 2, 3, and 4. Again, by defining the time where the correlations exceed a certain value C thres , we can define contours that indicate light-cones (see figure 4 with C thres = 0.01, 0.025, and 0.05). In general we observe a drastic change in the dynamics of correlations. While in the case of short-range interactions with α = 3, 4, we see a clear light-cone-like propagation behavior, for α = 1, 2, this behavior breaks down rapidly and instead almost instantaneous propagation of correlations is observed. This is summarized in figure 5a (left panel) where we show the light-cone boundary for C thres = 0.05. In the case of α = 3 and α = 4 we can extract a power-law exponent for the C thres = 0.05 contour, t C thres ∝ j η , with η = 0.67 ± 0.06 and η = 0.89 ± 0.02, respectively (dotted lines in figure 5a left panel).
Interestingly, the rapid change in propagation behavior is also directly reflected in the time-evolution of the experimentally much more accessible observableŜ x . Linked to the disappearance of the light-cone (figure 5 left panel) at α = 2, the behavior of the contrast decay as a function of time (figure 5 right panel) changes. For α = 3, 4 after an initial quadratic decay, Ŝ x decays slowly (surprisingly more slowly for α = 3 than for α = 4). For α ≤ 2, however, these two time-scales disappear and Ŝ x  exhibits a qualitatively different decay. The behavior for α ≤ 2 can be understood semi-analytically. For such long-range interactions, we can approximately replace the couplings by a constant: J ij ≈ J eff . Then, the XY Hamiltonian becomes a collective spin-Hamiltonian ∝Ŝ 2 x +Ŝ 2 y , whose dynamics is the same as that of the Ising Hamiltonian −J effŜ 2 z due to the conservation of total spin S 2 . The effective coupling constant J eff in our finite square lattice (M spins in total) can be determined, for example, by requiring that the total energy of the central spin interacting with all other spins with couplings J ij is the same as with coupling J eff . From the solution of the Ising model, one thus obtains S x (t)/M = cos M −1 (2tJ eff ) which is shown as dashed lines in figure 5a on the right. Given that there is some arbitrariness in defining a precise value for J eff , we see that the contrast decay for long-range interactions is fairly captured by this simple model.
An important question to answer is whether we can trust our dTWA predictions that indicate a critical behavior at α = D, i.e. when the power-law exponent matches the system dimension. This sharp transition is justified by the fact that for α < D the interaction energy becomes overextensive (diverges in the thermodynamic limit) and no bounds on the line-cone propagation can be established [32][33][34][35][36][37][38]. We test these predictions with exact solutions in a smaller 4 × 5 system (figure 5b) as well as for a 1D 21 × 1 chain (figure 5c). In the 1D case we calculate the correlations from the system center r i,c = (11, 1). While in the small 2D system, the crossover in the correlation propagation behavior is barely seen, and also the qualitative change in Ŝ x is not significant, in the 1D case (which features larger correlation distances than the small 2D system) a similar qualitative change in Ŝ x starts to emerge at α = 1. Note that also here the α = 1.5 case decays slower than the α = 2 case. In figures 5b,c we consider both, exact (points) and the dTWA solution (lines). Both, the exact and the dTWA solution feature the qualitative changes. Given the two facts: i) that in the exact 1D system calculations (with moderate size) a qualitative similar behavior to the one seen in the large 31 × 31 system starts to emerge ; and ii) that the agreement between the exact results and the dTWA solution becomes better in 2D, we are confident that the transitional behavior at α = 2 seen in figures 4 and 5a is indeed a real effect. Since the latter is also reflected in the decay of the Ramsey contrast, there is a good chance that future experiments will be able to verify the dTWA predictions in this regime.

Conclusion & Outlook
We have used a new numerical technique, the dTWA, to study the propagation of correlations in large 2D XY spin models with long-range interactions, in regimes accessible to current state-of-the art experiments with polar molecules or trapped ions. We benchmarked this new method in exactly solvable limits (Ising interactions and small systems) and found excellent agreement. In large systems, our method predicts a sharp change in the two-point correlations and contrast dynamics when the decay exponent crosses α = 2. Since other numerical techniques become inefficient in this regime, ultimately a quantum simulation experiment should verify these predictions.
In the future it will be interesting to study the nature of this sharp transition, not only in 2D, but also in 3D. Actually in 3D it is where polar molecule experiments using rotational states and exhibiting 1/r 3 interactions could observe this sharp change in the speed propagation of correlations. In this case it will be intriguing to investigate the role played by the anisotropic character of the interactions and the finite filling fraction on the propagation velocity of correlations. Systems where retardation effects in the dipolar interaction become relevant and 1/r 2 and 1/r terms can not be neglected (for example with atoms in two electronic states [20,58] or molecules prepared in vibrations states) could also become excellent laboratories for the observation of the dTWA predictions. In many realistic implementations, dissipation effects, induced for example by single photon spontaneous emission or cooperative radiation, can compete with the Hamiltonian dynamics. For modeling these experimentally relevant situations it will be important to adapt our technique to deal with a master equation formulation instead of pure Hamiltonian evolution.