Many-body dynamics of dipolar molecules in an optical lattice

Understanding the many-body dynamics of isolated quantum systems is one of the central challenges in modern physics. To this end, the direct experimental realization of strongly correlated quantum systems allows one to gain insights into the emergence of complex phenomena. Such insights enable the development of theoretical tools that broaden our understanding. Here, we theoretically model and experimentally probe with Ramsey spectroscopy the quantum dynamics of disordered, dipolar-interacting, ultracold molecules in a partially filled optical lattice. We report the capability to control the dipolar interaction strength, and we demonstrate that the many-body dynamics extends well beyond a nearest-neighbor or mean-field picture, and cannot be quantitatively described using previously available theoretical tools. We develop a novel cluster expansion technique and demonstrate that our theoretical method accurately captures the measured dependence of the spin dynamics on molecule number and on the dipolar interaction strength. In the spirit of quantum simulation, this agreement simultaneously benchmarks the new theoretical method and verifies our microscopic understanding of the experiment. Our findings pave the way for numerous applications in quantum information science, metrology, and condensed matter physics.

Advances in trapping and cooling ground-state polar molecules have produced ultracold, nearly degenerate gases [1][2][3][4], allowed complete control of their hyperfine, rotational, vibrational, and electronic degrees of freedom [5,6], and enabled their preparation in optical lattices [7]. Long-range dipolar interactions among molecules facilitate the exploration of fascinating manybody phenomena [8][9][10][11][12] and are useful for quantum information processing [13]. By encoding a spin-1/2 degree of freedom in two rotational states, Ref. [14] reported the first observation of dipolar-exchange interactions between molecules in an optical lattice, with signatures such as density-dependence of the spin coherence dynamics as probed by Ramsey spectroscopy (Fig. 1). Although such dynamics is expected to be governed by a particular spin model [15,16], a quantitative demonstration that the experimental observations are consistent with this model was lacking. A key obstacle was the challenge of modelling three-dimensional, strongly correlated, farfrom-equilibrium systems with long-range interactions.
In this work, we demonstrate that the molecules' Ramsey dynamics are quantitatively described by an XY spin model with anisotropic 1/r 3 interactions and rule out multiple alternative models. This validation required parallel experimental and theoretical developments. Experimentally, we demonstrate that we can tune the dipolar interactions (by a factor of two) by using two different rotational state pairs, and we systematically measure the dynamics under controlled conditions to quantitatively extract dependences on molecule number. Theoretically, we develop a new cluster expansion technique, which we term the "moving-average cluster expansion" (MACE), illustrated in Fig. 1(d), that is capable of describing the relevant non-equilibrium dynamics of spin models with long range and anisotropic interactions. MACE converges faster than the prior state-of-the-art disjoint cluster expansion (DCE) and thus enables the description of the spin dynamics of large numbers of molecules. We benchmark the convergence of MACE in a related (Ising) by encoding the spin in two rotational states [15][16][17].
(Alternative proposals to implement spin models in polar molecules are discussed in Refs. [18][19][20][21][22][23][24][25][26].) Here, S ± i and S z i are spin-1/2 operators satisfying [S z i , S ± j ] = ±δ ij S ± i , and the sums run over all occupied lattice sites with positions r i in units of the lattice spacing a = 532 nm. The dipolar interaction couples spins i and j with strength where Θ ij is the angle between r i − r j and the quantization axis, which in our system makes an angle of 45 • with theX andŶ lattice directions [see Fig. 1(b)]. The "exchange" or "XY" terms, S + i S − j + h.c., swap the spin states of molecules i and j. These arise from the transition dipole between the |↑ and |↓ rotational states, and they allow one molecule to flip from up to down while the other flips from down to up, together conserving their combined rotational energy. The "direct" or "Ising" terms, S z i S z j , arise because generally (at finite field) the dipole moments for |↑ and |↓ differ, and thus the parallel and antiparallel config-urations of two molecules have different energies. Note that Eq. (1) does not assume unit filling or even homogeneity; we will consider the details of the experimental distributions of molecules later. Figure 1 illustrates our experimental system, similar to that described in Ref. [14] (Methods and Supplementary Material contain details). We create N = 6, 000 to 23, 000 ground-state fermionic molecules of 40 K 87 Rb in the lattice. These molecules have lifetimes exceeding 25 s in the deep 3D optical lattice [7], with the lifetime limited by off-resonant light scattering. In our experiment we work at zero electric field, where the dipole moments and thus J z vanish, while the resonant exchange coupling J ⊥ remains finite. However, in order to benchmark the cluster expansions we will also theoretically study the case J ⊥ = 0, J z = 0.
To probe the spin system described by Eq. (1), we use the Ramsey protocol illustrated in Fig. 1(c). The spins are uniformly rotated to an equal superposition state by a resonant π/2 microwave pulse about theŷ spin axis, and at a later time t we read out the evolution of the spins by application of a final π/2 pulse that is phase-shifted by ϕ. This final pulse is equivalent to rotation about a vectorn = (sin ϕ, cos ϕ, 0) and measures the quantity cos ϕ i S x i −sin ϕ i S y i . By varying ϕ, we determine the global Ramsey fringe contrast C defined here as We include a π spin-echo pulse aroundŷ at time t/2 to remove the dephasing associated exclusively with inhomogeneous light shifts and isolate the effects of spin-spin interactions.
We vary the strength of the dipolar interactions by choosing different pairs of rotational states to realize the spin-1/2 system. The exchange coupling J ⊥ is determined by the transition matrix dipole element d ↑↓ = ↓ |d| ↑ , with d the appropriate spherical component of the dipole operator, via J ⊥ = −d 2 ↓↑ /4π 0 a 3 where 0 is the free space permittivity [16]. For the states used in Ref. [14], |↓ = |N = 0, m N = 0 and |↑ = |1, −1 , and for our lattice geometry, the transition dipole moment is expected to give an exchange frequency of |J ⊥ /(2h)| ≈ 52 Hz, where h is Planck's constant. For the alternative choice of rotational spin states, |↓ = |0, 0 and |↑ = |1, 0 , |J ⊥ /(2h)| is predicted to be twice as strong (∼100 Hz). This enhancement occurs because the molecules are aligned and oscillate along the quantization axis, as opposed to rotating about it, so that the dipole coupling is not reduced through time-averaging [16]. Figure 2 compares the observed contrast dynamics as a function of evolution time t for three different molecule numbers, N ∼ 5 × 10 3 , 1 × 10 4 , and 2 × 10 4 , and both pairs of rotational states. A larger coupling is observed for the {|0, 0 , |1, 0 } pair, which is promising for future applications of cold dipolar molecules. In particular, the faster coherent spin dynamics reduces the effects of technical limitations on the coherence times. The ability to compare theory and experiment for different exchange couplings is a key feature in our validation of the quantum simulation of the spin model in Eq. (1).
Theory.-To simulate the dynamics, we employ a novel cluster expansion technique. We compare it with a prior cluster expansion, the DCE, which breaks the system of spins into non-overlapping (disjoint) clusters, as shown schematically in Fig. 1(d). The algorithm to choose these clusters is given in Refs. [27,28]; roughly, it is designed to minimize the inter-cluster couplings. Observables such as the collective spin vector within each cluster are exactly computed numerically and then summed over all clusters to obtain the total collective spin dynamics.
The new approach developed here is to compute S α i (t) (α = {x, y, z}) by building an optimal cluster for spin i containing the spins connected to it by the g largest coupling constants V ij -see Fig. 1(d). The dynamics of S α i (t) is then computed by exactly solving the dynamics of the entire associated cluster by numerical integration, and the contrast is obtained by summing these spin expectation values over all i. We term our method the moving average cluster expansion (MACE). A natural advantage of the MACE method is that it is more robust to artifacts arising from finite cluster sizes; by constructing an optimal cluster for a each spin, MACE reduces surface effects wherein the dynamics of boundary spins in a DCE cluster may not be accurately captured.
We choose a molecule distribution according to our rough expectations of the experimental distribution based on the molecule formation process. In particular, the molecules are produced only at sites of the lattice initially populated by exactly one Rb and one K atom. Guided by the initial atomic numbers, temperatures, and trapping parameters (see Supplementary Information), we expect a doubly occupied Mott insulator domain of Rb atoms in the centre of the trap, surrounded by a unit-filled ellipsoidal Mott shell where molecules may be formed. Therefore, we assume a shell of molecules with a filling probability on site i given by , with α = 7 reflecting the ratio of axial and radial trapping frequencies, and the central shell radius of R c = 35. We determine the filling fraction by varying the shell width w, and set it to match the experimentally observed spin-echo contrast decay time for N = 1.2×10 4 molecules. Using this procedure, we find w = 30. Although our calculations use this specific distribution, we find that for a fixed local peak filling of molecules our theoretical results are largely independent of the chosen geometry.
Comparing theory and experiment.- Figure 2 demonstrates that our calculations (solid lines) quantitatively agree with the measurements for both rotational state choices and for a broad range of densities and evolution times. We emphasize that our results use a single global parameter to reproduce all the experimental data. This parameter is the ratio of the molecule filling factor f to the molecule number N .
We note a few interesting trends that emerge in both the theory and experimental data, then describe in the following paragraphs each of these trends in more detail. First, as in Ref. [14], we observe clear oscillations of the contrast. These oscillations are roughly independent of the molecule number N , but the frequency is found to be larger for the |1, 0 data, consistent with the enhanced spin-exchange coupling. Second, we observe that the spin coherence time decreases with an increase in lattice filling. This is a clear signature of spin-spin interactions. In comparing panels (a-c) with (d-f) of Fig. 2, this coherence time is also seen to be shorter for the spin states with larger spin-exchange coupling (|↑ ≡ |1, 0 ). Lastly, we find a trend of increasing oscillation amplitude for increasing molecule number.
The contrast oscillations arise from dipole-dipole interactions between the molecules. Clear oscillations at frequency |J ⊥ /(2h)| are visible for all measurements and calculations. Additionally, the Fourier transform of the theoretical contrast (Fig. 2, inset) clearly shows multiple oscillation frequencies on a broad, structured background. These frequencies are roughly determined by the size of the strongest couplings. The most prominent contributions apparent in the time evolution data appear at the frequencies ν/ √ 2, and ν/2 (next nearest neighbor coupling strengths), where ν = 104 Hz for |↑ = |1, 0 data and 52 Hz for |↑ = |1, −1 . In the Supplementary Materials, we have analyzed the experimental data by fitting to functional forms that oscillate at either a single frequency ν or at three frequencies. Based on a global analysis of fifteen data sets, we find clear evidence that a multi-frequency fit better captures the observed dynamics. Moreover, the analysis suggests a fundamental frequency of ν ∼ 108 Hz for the |↑ = |1, 0 data (reduced by half for the |1, −1 ), in excellent agreement with the expected value. The influence of the choice of rotational states is further illustrated in Fig. 3(a). This plot overlays contrast dynamics for the two pairs of rotational states, with the times of the |↑ = |1, −1 data rescaled by a factor of 1/2 to account for the different dipolar interaction strength. The collapse of the two datasets highlights that all of the observed dynamics arise from the dipolar interactions.
To investigate the coherence time τ 's dependence on the particle number N , we experimentally reduce the number of molecules using single-particle loss due to offresonant light scattering. This leaves the distribution of molecules invariant by reducing the density of particles uniformly (see Methods). We extract coherence times by fitting the contrast dynamics to exp (−t/τ ) for both the theory and experiment; the results are shown in Fig. 3(b). Whereas the oscillations come mainly from the largest frequencies (smallest spacings) between molecules, the decoherence arises from interactions of molecules at a variety of spacings. We expect τ ∝ 1/N , which is a characteristic signature of interactions. The scaling arises because, as the lattice filling f ∝ N increases, the mean distance between molecules decreases asR ∼ f −1/3 , leading to an average dipolar interaction that scales as 1/R 3 ∼ f . We indeed see an approximate scaling τ ∝ 1/N , as well as a factor of 2 reduction in coherence time for the |1, 0 data compared to the |1, −1 data. We note that the theory here uses a peak filling of f = 8% for a molecule number of N = 2×10 4 , which is within a factor of two of estimates based on loss measurements (∼ 9% for N ∼ 1 × 10 4 ) and direct imaging [14,29].
Finally, the oscillation amplitude increases with N for the experimental fillings studied since the probability of a molecule having an occupied nearest neighbour site increases with N . To characterize the amplitude of the oscillations, we plot in Fig. 3(c) the root-mean-square (rms) residuals of the data from the exponential fit. We find that for increasing molecule number N there is a systematic increase of the residuals, due to oscillatory dynamics absent in the simple exponential fit.
Discussion.-Given the accuracy of our theoretical cluster expansion (see Methods and Supplementary Information), our measurements can rule out multiple alternatives to the spin-exchange model Eq. (1) for describing the experimental observations. For example, the experimentally measured dynamics is inconsistent with the Ising model, where the contrast oscillation amplitudes are significantly smaller than those observed experimentally (see Methods and Fig. S2 of the Supplementary Information for a demonstration.) As another example, Fig. 4(a) shows that the spin-exchange model with interactions truncated to the nearest neighbours in the lattice also fails to reproduce our measurements even qualitatively.
We can now make predictions for future experiments, for example the effects of increasing the filling and varying the initial spin preparation. Figure 4(b) shows the effects of increasing the filling towards unity. The contrast decay becomes even more rapid. The contrast oscillation amplitude initially increases monotonically with increasing filling (see Fig. 4) and then saturates around fillings of 25%, and finally decreases. At even higher fillings, the dynamics becomes overdamped and the contrast oscillation amplitude decreases, similar to the one-dimensional behavior studied in Ref. [12].
We find an interesting feature when we examine the dependence on the initial spin tipping angle. For arbitrary tipping angles θ = π/2, the rescaled contrast C(t)/ sin θ is largely independent of θ, as shown in Fig. 4(c). The contrast decay arises from two effects: dephasing of precession angles between spins, since the mean-field precession rate differs spin to spin (this only plays a role for θ = π/2), and contrast decay due to the buildup of correlations. If ij V ij = 0, as it nearly does for the anisotropic dipolar interaction, one can show using the analytic formulas for short-time dynamics in Ref. [12] that these two terms balance in such a way that C(t)/ sin(θ) is independent of θ. It is very interesting that the collapse of the curves in Fig. 4(c) holds to longer times. Our theory also predicts that the phase φ = arcsin [ i S y i / i S x i ] diffuses randomly and remains small, rather than undergoing a mean-field-like precession, as shown in Fig. 4(d).
Summary and conclusions.-We have experimentally measured the Ramsey contrast dynamics in ultracold KRb molecules, varying both the molecule density and the strength of the dipolar interactions. We have developed a theoretical method that allows us to efficiently and accurately compute the spin coherence evolution of macroscopic numbers of molecules. The theory quantitatively reproduces the experimental data, and our comparison allows us to determine a typical experimental lattice filling fraction between 5-10%, which is roughly consistent with that obtained by independent measurements [29]. One question for future investigation is the role of inhomogeneous light shifts, which add a spatially varying i ∆ i S z i term in Eq. (1) and can suppress exchange interactions. If we use values of ∆ i expected from AC polarizabilities [6], our data and theory begin to disagree at times longer than about 20 ms. This could arise from mis-estimates of the polarizabilities, spatial location of the molecules, or from unanticipated technical sources of decoherence.
The results presented here exemplify how a clean ultracold system can guide the development of theories applicable to the numerous quantum systems governed by similar Hamiltonians. For example, MACE can help understand recent experiments using room-temperature polycrystalline molecular solids [30], or atomic optical lattice clocks [31,32]. We expect the MACE method to have applications to numerous other systems described by long-range spin models [33][34][35][36][37][38][39][40][41], such as magnetic atoms, molecules, Rydberg atoms, trapped ions, solidstate defects, as well as exciton transport in molecules, relevant for chemistry and biology [42,43].
In the future it will be fascinating to examine the development of correlations more directly, both theoretically and experimentally, and thereby explore transport/thermalization or lack thereof, e.g. glassiness and many-body localization [44][45][46][47]. Acknowledgements

Experimental production of ultracold KRb in an optical lattice
To produce the molecules in a lattice, we follow the experimental procedure detailed in Ref. [14]. Briefly summarizing, we create between N = 6, 000 to 23, 000 ground-state molecules in a lattice of spacing a = 532 nm. We start from 40 K atoms initially at T /T F ≈ 0.5 and 87 Rb atoms initially at T /T c ≈ 0.5 (T F and T c are the Fermi and Bose-Einstein condensate transition tempera-tures, respectively), confined in traps with 25 Hz radial and 185 Hz axial frequencies for Rb. We create Feshbach molecules by ramping the magnetic field across a broad interspecies Feshbach resonance at ∼546 G, and then transfer these to ground-state molecules by a stimulated Raman adiabatic passage (STIRAP) process. To realize spin systems, we freeze out the tunnelling by working in a deep 40 E r lattice, where E r =h 2 k 2 /(2m) is the photon recoil energy with momentum k = π/a, where m is the molecule mass. We image by reversing the STI-RAP process and then taking absorption images of the K atoms after release from the lattice.
In order to change the filling of molecules uniformly, leaving their distribution unaltered, we employ single particle loss from light scattering. Specifically, we hold the molecules in the deep lattice with two additional running beams that are homogeneous over the scale of the cloud. The time for which these are applied is chosen to achieve the desired density.

Convergence of theory
We demonstrate the convergence of our theoretical method (MACE), validating our theoretical calculations. We also discuss the advantages of the MACE over prior methods (DCE). We demonstrate convergence for the spin-exchange dynamics, but to further verify the accuracy and convergence to the correct solution, we also consider the J ⊥ = 0 case, which admits an exact analytic treatment [12,[48][49][50]. For example, for an initial tipping angle of the spins relative to −ẑ, θ = π/2, the exact Ramsey contrast is C Ising = i j =i cos (J z V ij t/2) (h = 1). Figure S2 in the Supplementary Information illustrates the convergence properties of the DCE and MACE. Dashed (solid) lines show the contrast versus time for the DCE (MACE) with cluster sizes g = 4, 6, 8, 10, in both cases for N = 5000 spins in a ellipsoidal shell with prob- with R c = 25, α = 3, and w = 25. Figure S2(a) shows the computed XY dynamics (J z = 0), while the panel (b) shows Ising dynamics (J ⊥ = 0) compared with the exact Ising solution (open circles). In both cases, the MACE is quantitatively accurate even for small clusters g ∼ 6. In contrast, although the oscillations at shorter times are well-captured by the DCE, the long time behaviour is poorly converged, even for g = 10. For the Ising case, despite the apparent convergence of the DCE we find that a large discrepancy remains even for enormous cluster sizes of g ∼ 500. These statements are complemented in the Supplementary Information by a quantitative analysis of the root-mean-square error, both between successive cluster sizes and for the Ising case between the approximate and exact solution. The reason the DCE fails is that the dynamics of spins at cluster boundaries is captured poorly, while in the MACE each cluster is optimized around the molecule contributing to the spin expectation, avoiding these problems. Although for sufficiently large clusters these boundary contributions are negligible and both methods formally converge to the same result, the cluster sizes required by the DCE to reach this limit are far larger than numerically possible. We also note that the MACE is trivially parallelizable: each spin can be considered independently. In contrast, the creation of clusters in the DCE is less trivial to parallelize and can hinder scaling to large numbers of spins.

DIFFERENTIAL LIGHT SHIFT IN A 3D LATTICE
To enable our use of coherent Ramsey spectroscopy on the sample of ultracold molecules, we create a nearly "magic" lattice trap for the two spin states |↑ and |↓ , i.e. one that is nearly spin-state-independent [14]. Because the AC polarizability of molecules is anisotropic [6], the angle between the polarization of the light fields and the quantization axis can be used to tune the light shifts for different molecular states. In Ref. [14], the details of reducing the differential AC polarizability for the case of |↑ ≡ |1, −1 were reported. The polarizations of thê X,Ŷ , andẐ lattice beams were chosen to be linear, at angles of +45, −45, and +45 degrees with respect to the quantization axis, respectively. In this work, we additionally study the case |↑ ≡ |1, 0 . In this case, the smallest differential AC polarizability is achieved by setting the linear polarizations to the "magic" angle for the states |0, 0 and |1, 0 , as determined in Ref. [6]. We find the differential light shift to be minimized when the linear polarizations of the three lattice beams are at angles of roughly +54, −54, and +54 degrees with respect to the quantization axis, respectively.

EVIDENCE FOR MULTIPLE INTERACTION ENERGIES
A key feature of the observed spin coherence dynamics is oscillations due to dipole-dipole interactions. While previously spin coherence data has been fit to a function with a single oscillation frequency [14], many interaction frequencies are expected due to the long-ranged and anisotropic nature of the dipolar interactions. Here we more comprehensively compare the data to empirical functions that include either a single frequency or multiple frequencies. We simultaneously fit over a dozen data sets that include both sets of spin states, i.e. |↑ = |1, −1 or |↑ = |1, 0 , rescaling the time axis for the |1, −1 data by a factor of 1/2 to account for the known difference in dipolar couplings. In considering a single frequency, we fit to the functional form A cos 2 (νπt)+(1−A) exp(−t/τ ), which fixes the contrast to unity at time t = 0. In considering multiple frequencies, we use our knowledge of the molecules' geometry in the lattice and of the form of the interaction to set the frequency ratios to match those of the 3 strongest pairwise interactions, and fit to the form i A i cos 2 (ν i πt) Fig. S1 (a), we show the results of a single-frequency and three-frequency fit for a single data set with |↑ = |1, 0 . The fit for multiple frequencies clearly does a better job at capturing the observed dynamics, as compared to the single-frequency fit. To make this statement more quantitative, taking into account the additional free parameters allowed for by the multi-frequency fitting, we analyze the cumulative reduced chi-squared (χ 2 r ) of the fits as applied to all the data sets. The result of this analysis is shown in Fig. S1 (b), as a function of the assumed oscillation frequency. Two major pieces of information can be gleaned from this analysis. First, the fitting with multiple frequencies does a much better job of describing the observed contrast oscillation data, and yields a considerably lower reduced chi-squared value. Second, we find from the frequency dependence of χ 2 r that the observed oscillations are most consistent with a primary oscillation frequency of ∼ 108 Hz for the |1, 0 data (and by construction a factor of two lower for the |1, −1 case). This value is in good agreement with the theoretical prediction of 104 Hz. Figure S2 demonstrates the superior convergence of the MACE method, and is discussed in detail in the Methods section of the main text. We also study an alternative quantitative measure of the two cluster expansions' convergence: we show (1) for both XY and Ising, the root-mean-square difference between adjacent cluster sizes and (2) for Ising, the root-mean-square error to the exact solution. That is, for both cases, we calculate ∆(g) = (1/T ) T 0 dt (C g+1 (t) − C g (t)) 2 where C g (t) is the contrast at time t calculated for cluster size g. For the error to the exact solution, the same formula is used with C g+1 being replaced by the exact solution.

BENCHMARKING THEORETICAL CONVERGENCE AND ACCURACY
For the current comparison, we use an integration time T = 30h/J, with J = J ⊥ or J z as appropriate. The ∆(g) versus g are shown in log-log plots in Fig. S3. Open symbols are adjacent differences, while closed symbols are the difference to the exact theory. The dramatically lower errors and faster convergence of the MACE are evident.