Synthetic quantum systems with interacting constituents play an important role in quantum information processing and in explaining fundamental phenomena in many-body physics. Following impressive advances in cooling and trapping techniques, ensembles of ultracold polar molecules have emerged as a promising platform that combines several advantageous properties [1][2][3][4][5][6][7][8][9][10][11] . These include a large set of internal states with long coherence times [12][13][14][15][16][17] and long-range, anisotropic interactions. These features could enable the exploration of intriguing phases of correlated quantum matter, such as topological superfluids 18 , quantum spin liquids 19 , fractional Chern insulators 20 and quantum magnets 21,22 . Probing correlations in these phases is crucial to understanding their properties, necessitating the development of new experimental techniques. Here we use quantum gas microscopy 23 to measure the site-resolved dynamics of quantum correlations of polar 23 Na 87 Rb molecules confined in a two-dimensional optical lattice. By using two rotational states of the molecules, we realize a spin-1/2 system with dipolar interactions between particles, producing a quantum spin-exchange model 21,22,24,25 . We study the evolution of correlations during the thermalization process of an out-of-equilibrium spin system for both spatially isotropic and anisotropic interactions. Furthermore, we examine the correlation dynamics of a spin-anisotropic Heisenberg model engineered from the native spin-exchange model by using periodic microwave pulses [26][27][28] . These experiments push the frontier of probing and controlling interacting systems of ultracold molecules, with prospects for exploring new regimes of quantum matter and characterizing entangled states that are useful for quantum computation 29,30 and metrology 31 .

Ultracold dipolar particles trapped in optical arrays have emerged as a powerful and flexible platform to probe idealized models of quantum many-body physics with long-range interactions 32,33 . For example, studies that used polar molecules [9][10][11] , Rydberg atoms 34 and magnetic atoms 35 have explored a host of quantum phenomena including quantum magnetism 25,28,36 , symmetry-protected topological phases 37 and extended Hubbard models 38 . The platform of polar molecules in particular has several unique properties [9][10][11] . Ultracold polar molecules possess a tunable dipole moment in their electronic ground state, enabling strong, long-range interactions compatible with negligible population relaxation over experimental timescales. In addition, compared to atoms, molecules have rotational and vibrational degrees of freedom that can be harnessed to store information in a large local Hilbert space and to engineer a rich set of many-body Hamiltonians. However, this complex internal structure also presents its own experimental challenges in the preparation and detection of quantum-state-controlled ensembles of molecules.

Two complementary approaches have been developed to address these challenges. On the one hand, the assembly of heteronuclear molecules from evaporatively cooled atomic quantum gases 1 has generated bulk ensembles with high phase space density, enabling the recent creation of degenerate molecular Fermi gases for which the collective properties can be characterized with absorption imaging techniques 3,7 . On the other hand, molecules have been prepared in optical tweezers 39 , realizing a bottom-up approach to controlling small numbers of molecules with single-molecule detection capabilities 4,5 . Here we demonstrate quantum gas microscopy of polar 23 Na 87 Rb molecules in an optical lattice, combining features of both approaches: large numbers of molecules in their hyperfine, rovibrational and motional ground state, coupled with single-molecule detection. Quantum gas microscopy enables high-fidelity and simultaneous optical detection of ensembles of particles in a regime of low temperatures, high density and strong interactions. Pioneered with atomic gases, this technique has enabled unprecedented local observations of quantum phase transitions, spin and charge correlations in Hubbard systems and impurity physics 23 . Recently, we have extended microscopy techniques to quantum gases of non-interacting excited-state molecules, observing bosonic bunching correlations in their density fluctuations 40 .

In this work, we apply our ability to measure microscopic correlations to study out-of-equilibrium spin systems realized with interacting molecules. To that end, we prepare the molecules in their rovibrational ground state, which possesses a large body-frame electric dipole moment. We encode a spin-1/2 in the ground and first excited rotational states of the molecules. An effective spin interaction arises from the resonant exchange of rotational excitations between pairs of molecules, as was previously demonstrated with polar 40 K 87 Rb molecules 25,41 . Here, with molecules pinned in a deep two-dimensional (2D) optical lattice, the system realizes the quantum XY spin-exchange Hamiltonian 21,22,24 r r

are spin-1/2 operators for molecule i, r i is the position of molecule i in units of the lattice constant a lat and θ is the angle between the quantization axis and the vector a. The quantum average ⟨•⟩ accounts for the finite size of the molecule centre-of-mass wavefunctions. XY Hamiltonians are some of the most extensively studied models for magnetism and excitation transport. However, like most quantum spin models, many open questions in computationally intractable regimes remain, such as the prediction of low-temperature phases in the presence of frustration or the long-time dynamics of disordered or driven systems. Polar molecules are poised to be a leading platform to simulate such spin models because of the large ratio of spin coherence times to the interaction timescale h/| J| that we demonstrate in this work, which compares favourably to other systems used for studying similar physics such as Rydberg atoms 28 . We observe the effects of the dipolar interactions by performing a series of quench dynamics experiments, in which we initialize the molecules in a product state and then evolve under the interaction Hamiltonian in equation ( 1), after which we measure the correlations between the molecules with single-site resolution. For a lattice with low filling, the correlations exhibit weakly damped oscillations with a frequency that depends on the site displacement, directly reflecting the long-range and anisotropic character of the interactions.

We start our experiments by preparing an array of ground-state NaRb molecules 42 in a 2D optical lattice (Fig. 1a). The molecules are formed by loading degenerate Bose gases of Na and Rb into the lattice and using magnetoassociation followed by stimulated Raman adiabatic passage (STIRAP) to convert atom pairs into molecules in their rovibrational and hyperfine ground state 43,44 . The lattice is sufficiently deep such that tunnelling of molecules between lattice sites is negligible over the timescale of the experiments. The maximum achieved filling of ground-state molecules in the lattice is 15(1)%, taking into account the ground-state molecule detection efficiency of 92.0(3)% (Fig. 1c).

To simulate an effective spin-1/2 system, we use microwaves near 4.18 GHz to drive transitions between a hyperfine state in the ground manifold (N = 0) and a hyperfine state in the first excited rotational manifold (N = 1), labelled |↑⟩ and ↓⟩ | , respectively. For imaging, we transfer molecules in |↑⟩ to the weakly bound Feshbach state, which we then dissociated. We detect the resulting Rb atoms with fluorescence imaging, allowing us to extract the position of each molecule in |↑⟩ with single-site resolution 40 .

To benchmark the coherence of the two-level system encoded in the rotational states of the molecules, we perform a Ramsey interferometry experiment. Starting with all of the molecules in ↑⟩ | , we apply two π/2 pulses separated by a total free precession time T (Fig. 2a). A spin echo π pulse can be inserted in the middle of the precession time, which mitigates dephasing owing to quasi-static sources of single-particle to extract the contrast A/B. We dilute the sample to have a peak molecule filling of 1.0(2)% to reduce the influence of dipolar interactions between the molecules (additional data at higher densities are shown in Extended Data Fig. 3).

The results of this experiment are shown in Fig. 2. Without a spin echo pulse, the observed dephasing (1/e time of 56(2) ms) is faster than that predicted for an interacting XY spin system described by equation (1), most probably because of the residual differential light shifts (Methods). However, with the addition of a single spin echo π pulse, we observe that the decay of the experimental contrast is consistent with the dynamics expected from equation (1) out to 400 ms. This indicates that the coherence time of the two-level system is much longer than the expected millisecond-scale nearest-neighbour interaction time, which is a favourable regime for experiments studying coherent many-body dynamics, such as spin-squeezing or preparing rotational superfluids 31,45 .

Having established that the molecules in the lattice can be modelled as a closed quantum spin system out to long times, we probe the growth of their correlations due to dipolar interactions after a quench. For these experiments, we set the peak molecule filling to 5.4(4)%, which is deliberately chosen to be lower than our maximum filling shown in Fig. 1 to slow down the thermalization timescale compared to the nearest-neighbour interaction timescale. Under these conditions, we expect to observe coherent oscillations in the dynamics of various observables. To set the quantization axis for the molecules, we apply an external magnetic field pointing perpendicular to the 2D lattice plane so that the dipolar interactions between the molecules are isotropic. The quench is initiated with a global π/2 rotation to transfer the molecules from |↑⟩ to X

. The interacting molecules then evolve for a time T, with an odd number of spin echo pulses introduced to mitigate residual effective fields (Methods). Finally, we rotate the spins into the measurement basis with a final π/2 pulse, which is 180° out of phase with the initial pulse so that in the absence of any many-body interactions all of the molecules are transferred to ↓⟩ | and, therefore, appear dark to our detection scheme. Any mole- cules in ↑⟩ | are then detected and their positions in the lattice are recorded.

The results of this experiment, along with comparison to numerical simulations, can be seen in Fig. 3. Initially there are a negligible number of molecules in ↑⟩ | , but by 18 ms approximately 30% of the molecules are in |↑⟩. However, this overall loss of magnetization during the evo- lution of the system is not monotonic, as clear oscillations in the molecule number can be seen in Fig. 3a in good agreement with the theoretically expected frequency of |V(e)|/2h = 241 Hz, where e = (1, 0), corresponding to the interaction energy between spins on neighbouring sites. Although such oscillations have been previously observed 25 , our site-resolved measurements allow us to extract the lattice-averaged correlation function

↑ yielding a more detailed characterization of the dynamics of the system. Here, ⟨•⟩ denotes averaging over the quantum state of the many-body spin system as well as classical realizations of the dilute lattice filling, and N s is the number of sites in the region used to evaluate the correlator. The normalization factor is N ρ = (⟨ ⟩ ) undetected state |↓⟩, whereas interactions play an important role in the dynamics of pairs or clusters of nearby spins, so they can evolve into the detected state in a correlated way. For example, the quantum state of an isolated pair of molecules with displacement a at the end of a Ramsey sequence with evolution time T is

V T ħ cos( ( ) /4 ) ↓↓ up to a global phase factor. This two-molecule entangling-disentangling dynamic is responsible for an oscillatory behaviour of the correlations given by C T V T ( , ) sin ( ( ) /4 ) 2 a a ℏ ∼ in the limit of vanishing lattice filling and for short evolution times. This limiting case underlies the correlation oscillations in the data in Fig. 3c. It was observed that these oscillations damp slowly both in the data and in exact diagonalization calculations by using the dipolar XY model. Because the numerics are performed for a closed quantum system, we understand the damping to be due to corrections with the simple two-spin picture above. Although at short times pairs of nearby spins are the dominant contribution to the dynamics, at later times couplings to spins that are farther away play an increasingly important role.

Another defining characteristic of the dipole-dipole interaction is its spatial anisotropy. By performing the same quench dynamics experiment but with the magnetic field in the plane and tilted 9° with respect to the lattice axes, we can clearly see the effect of the anisotropic interaction potential in the dynamics 46 . The resulting correlations are shown in Fig. 4, alongside the full correlation matrices for the isotropic configuration (discussed in the last section), to highlight the change from isotropic to anisotropic correlations. For example, in contrast to the correlations shown for the isotropic case, the correlations along the x axis grow faster than the correlations along the y axis (where x and y are the lattice axes), as expected because the interaction energy is a factor of approximately 2 stronger along this axis. In addition, the interaction potential is zero for θ = 54.7°, which falls near the y = x diagonal, and, therefore, we expect correlations to be highly suppressed along that direction. We observe this in the experimental correlation data in Fig. 4h, in which it can be seen at 4 and 5 ms that the y = -x diagonal has strong correlations but there are negligible correlations for the y = x diagonal. Our correlation data demonstrate the ability to easily tune the molecular interaction anisotropy with the direction of the magnetic field, which will be useful in future studies of frustrated spin systems.

Although the XY model is the natural spin Hamiltonian implemented in our polar molecule system, one can go beyond the resonant spin-exchange term to realize more complex interactions, such as the XXZ or XYZ models, by using microwave Floquet engineering. The essence of this technique is to modulate the Hamiltonian in time at a rate faster than the interaction timescale, so that the time-averaged Hamiltonian is different from the original Hamiltonian [26][27][28] . Here we use rapid periodic π/2 rotations to realize an effective XXZ Hamiltonian. An additional π pulse in the centre of the cycle serves as a spin echo to cancel the dephasing due to any effective fields without affecting the many-body interactions. The total Floquet cycle can be seen in Fig. 5a.

The effect of each Floquet cycle is to evolve the molecules for an equal amount of time under H XY and H XZ , time averaging to

which is equivalent to the XXZ model with a permutation of the axis labels.

One important feature of the Hamiltonian in equation ( 2) is that the magnetization along X is conserved, whereas that along Y and Z are not, owing to the U(1) symmetry corresponding to rotations about the X axis. Therefore, a key test of engineering the correct effective Hamiltonian is to verify this spin conservation. We prepare three distinct initial states, corresponding to all of the spins pointing along X, Y or Z, and then subject the spins to the same Floquet Hamiltonian by using the pulse sequence shown in Fig. 5a. The density of the molecules was 8.1(8)%, and the quantization axis was perpendicular to the 2D lattice so that the interactions are isotropic in the plane. Finally, we measured the total number of spins that remained in the initial state for different evolution times.

The results of this experiment are shown in Fig. 5b. The magnetization of molecules initially prepared in the X + ⟩ | state remains approximately unchanged throughout the dynamics, whereas that of molecules prepared in the Y + ⟩ | and Z + ⟩ | states decays towards the unpolarized value of 1/2. In addition, the data is Y-Z symmetric, as expected from equation (2). Comparing the experimental data to exact diagonalization calculations by using an XYY model shows reasonable agreement, indicating that the time-averaged Hamiltonian is an appropriate descriptor of the observed dynamics.

Finally, we can measure the molecule correlations as a function of time to test whether the oscillation frequency changes for the Floquet dynamics compared to the original XY Hamiltonian. Starting with a peak molecule filling of 3.1(4)%, we prepare all of the spins in ↑⟩ | , quenched the system to evolve under the Floquet Hamiltonian and then measure the correlations between the molecules after each Floquet cycle. By comparing the data in Fig. 5c to the XY experiment in Fig. 3c, we can see that the oscillation frequency for the nearestneighbour correlations has been halved, as expected for the desired Floquet Hamiltonian. To conclude, we have demonstrated that microwave Floquet engineering of polar molecule systems enables the study of spin-anisotropic Heisenberg models, further enriching the quantum simulation toolbox in this system.

The experiments described above demonstrate the powerful capabilities of ultracold molecules to study dynamics in the context of quantum magnetism, and create new possibilities for studying quantum physics with ultracold molecules more broadly. The ability to detect correlations with single-site resolution could be used in other experiments to distinguish different equilibrium phases of matter, such as rotational superfluids that should emerge by adiabatic preparation of the ground state of the XY Hamiltonian even at the current demonstrated filling fractions 45 . The single-site resolution also enables optical manipulation of the underlying potential at the smallest relevant length scale, allowing the initialization of arbitrary configurations to study impurity physics or transport dynamics 23,47 . In addition, although so far we have demonstrated interactions between molecules that are frozen in the lattice, complex phases of matter are also expected to emerge upon allowing the molecules to tunnel throughout the lattice. The interplay between the dipolar interactions and the kinetic energy of the molecules is predicted to yield a rich phase diagram consisting of superfluid and supersolid phases, as well as fractional Mott insulators 24,48,49 . Our platform is particularly well-suited to study this physics owing to the strong dipolar interaction strength compared to magnetic atoms 35 , as well as the lack of collective decay mechanisms of the molecules compared to Rydberg-dressed atoms 46 . Finally, the capabilities demonstrated here to measure correlations between interacting polar molecules may prove useful in the context of quantum computing 29,30 , for which they could be used to verify the creation of Bell and Greenberger-Horne-Zeilinger states, as well as for quantum error correction.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05558-4.

The molecules are created by first loading a doubly degenerate mixture of Na and Rb atoms into a single plane of a three-dimensional (3D) optical lattice. The 2D lattice in the x-y plane has a 752 nm spacing and the vertical lattice has a 3.8 μm spacing. We then sweep the magnetic field across an interspecies Feshbach resonance at 347.6 G to form weakly bound molecules. Each lattice site that contains one Na and one Rb atom forms a weakly bound molecule 43,44 , whereas sites with more than one atom of each species are emptied by three-body loss processes. The remaining Na and Rb atoms that do not form molecules are removed with resonant light pulses. The 2D lattice depth for all of the experiments is 49E R with the exception of the coherence data for which the lattice depth is 34E R . Here E h ma

2 , where m is the molecule mass and a lat = 752 nm. At both depths, the tunnelling of the molecules in the lattice is negligible over the course of the experiment. Details of the optical potentials and Feshbach molecule creation process can be found in ref. 40 .

After forming the Feshbach molecules, we ramp the magnetic field to 335 G and transfer the molecules to the electronic and rovibrational ground state by using STIRAP. Following the procedure described in ref. 42 , we use two external-cavity diode lasers (Toptica DL Pro) with wavelengths of 770 nm and 1248 nm, locked to a common high-finesse ultra-low expansion cavity (Stable Laser Systems) by using the the Pound-Drever-Hall technique. The cavity free-spectral range is 1.5 GHz and the finesse is 34,000 at 770 nm and 43,000 at 1,248 nm. We use a fibre electro-optic modulator for each laser (EOSpace) to tune the frequency of the laser relative to that of a fixed cavity mode, as well as to generate sidebands for the Pound-Drever-Hall lock. The laser linewidth is narrowed by active feedback of the diode current, whereas slow drifts in the frequency are compensated with a separate feedback loop stabilizing the laser piezo. Each laser uses a FALC 110 (Toptica) laser locking module for both feedback loops. The Rabi frequencies are 2π × 0.70(4) MHz and 2π × 0.9(1) MHz for the 1,248 nm and 770 nm transitions respectively, and we typically achieve 93.9 . To deliberately dilute the density of the molecules for some of the experiments presented in the paper, we reduce the one-way STIRAP efficiency by decreasing the optical power in the 770 nm beam to reach the desired filling fraction.

To detect the molecules, we reverse the STIRAP process described above to convert molecules in |↑⟩ back to weakly bound Feshbach mol- ecules. The Feshbach molecules are dissociated and we then image the corresponding Rb atoms by using fluorescence to tag the positions of the molecules in ↑⟩ | within the lattice 40 . We do not detect the molecules in |↓⟩, so our current measurements do not distinguish between that state and an empty site, although this may be achieved in future work with bilayer techniques used for spin-resolved imaging in atomic microscopes 50,51 .

To simulate an effective spin-1/2 system, it is necessary to couple two states that are well separated in energy from other molecular states and that are as insensitive as possible to noise from variations in magnetic or optical fields. To achieve this condition, we use microwaves resonant with the transition between the ground and first excited rotational state of the molecules. In general the molecule can be in a superposition of states labelled N m m m , , , ⟩

, where N is the rotational angular momentum of the molecule, m N is its projection onto the quantization axis set by the external magnetic field, and m I,Na(Rb) is the projection of the nuclear spin of the Na(Rb) atom. A representative microwave spectrum with several of the rotational and hyperfine transitions is shown in Extended Data Fig. 1, along with a representative Rabi oscillation between the specific pair of states used in this work.

One important source of single-particle decoherence is the differential dynamic polarizability 13,16,[52][53][54][55] ). This causes an unwanted shift of microwave transition frequencies across the trap because of the spatially varying intensity profile of the trapping light. To mitigate this, we choose two specific magnetic field configurations that produce near 'magic' trapping conditions, where

For the experiments with isotropic interactions, the magnetic field is 60 G, whereas for anisotropic interactions the magnetic field is 4.1 G. Although a higher field would have been preferable for the latter case to increase the energy splittings to other states outside of the Hilbert space of the spin system, our coil geometry limits the maximum magnetic field that can be applied in the lattice plane. Our technique makes use of the fact that, at zero electric field and weak magnetic fields, m N is not a good quantum number for the molecular Hamiltonian owing to strong hyperfine couplings, especially the nuclear quadrupole moment coupling in N = 1. This coupling gives a specific N = 1 state an admixture of other hyperfine states, which can be leveraged to match the ↓⟩ | polarizability to the |↑⟩ polarizability. 54,56 . Calculated differential polarizabilities are shown in Extended Data Fig. 2. Over the spatial extent of the cloud, the microwave transition frequency varies by less than 20 Hz because of the inhomogeneous intensity profile. This is consistent with the observed coherence decay time in the absence of a spin echo pulse (Fig. 2).

After transferring the molecules to the rovibrational ground state, transitions to the first excited rotational state are induced using microwaves near 4.18 GHz. The microwaves are generated by mixing a local oscillator at 4.13 GHz provided by an analogue signal generator (Agilent E8257C) with an intermediate frequency provided by an arbitrary waveform generator (Keysight 33600A). The intermediate frequency is 50.13747 MHz when the magnetic field is 60 G, and 50.08300 MHz when the magnetic field is 4.1 G. The mixed signal is then amplified (Mini-Circuits ZHL-5W-63-S+) before being sent to a home built helical antenna mounted underneath the vacuum chamber. We programme the Keysight arbitrary waveform generator to set the duration, amplitude and phase necessary for each pulse within the experiments, and keep the Agilent signal generator at fixed frequency, amplitude and phase. The oven-controlled crystal oscillator inside the Agilent signal generator is used as a common 10 MHz clock for both instruments.

We use a microwave Rabi frequency of 9 kHz for all of the experiments with spatially isotropic interactions. This frequency is chosen to be large compared to the intermolecular interactions (<1 kHz) but small compared to the splitting between the different hyperfine states to avoid driving transitions outside of the effective two-level system. However, for the spatially anisotropic correlation measurement we decrease the Rabi frequency to 4 kHz to minimize Fourier broadening from the microwave pulses, since the hyperfine states are more closely spaced at the lower magnetic field that we use for these experiments.

In the Ramsey spectroscopy experiments shown in Fig. 2, we observe a slow drift in the phase of the fringe in addition to a decay in the fringe amplitude. This phase drift is not expected from the desired Hamiltonian in equation ( 1), and can be caused by a combination of the interactions and the inhomogeneous light shifts, or an unknown time-varying field present in the lab. We observe that adding multiple π pulses does not affect the decay rate of the fringe amplitude, but it does remove the phase drift. Therefore, for all of the dynamics experiments in Figs. 3 and4, we use one π pulse for any data collected with less than 10 ms of evolution time, for which the phase shift is negligible. For data collected after 10 ms of evolution time, more than one π pulse is required to mitigate the phase shift, so we use three π pulses during the evolution time, deliberately keeping the number of π pulses odd. In addition, for each experiment we synchronize the start of the first microwave pulse with the alternating current line to ensure that every experiment is performed with the same background magnetic field conditions.

The pulse sequence for the Floquet engineering protocol is shown in Fig. 5a. For these experiments, we set the Floquet evolution time to τ = 1 ms and use a π-pulse duration of 56.2 μs so that each Floquet cycle is short compared to the nearest-neighbour interaction period of the original XY Hamiltonian 2h/|V(e)| = 4.15 ms seen in Fig. 3, where e = (1, 0). In addition, for the experiments in which we measure the decay of the magnetization for the initial states | Y + ⟩ and | Z + ⟩ (Fig. 5b), an extra π pulse is added at the end of the evolution time to return the coordinate system of the Bloch sphere to its initial orientation before measurement.

In Fig. 2, we report a Ramsey spectroscopy measurement of the coherence of our two-level system at the lowest achievable lattice filling of 1%. Here we include the results from additional Ramsey spectroscopy experiments, which are performed at higher lattice fillings. This increases the influence of dipolar interactions in the system, allowing us to test the effect of these interactions on the coherence time. The experimental parameters are the same as that reported for Fig. 2, and a spin echo pulse is used to eliminate decoherence from quasi-static sources to focus primarily on the role of the dipolar interactions. We show the results for several lattice fillings in Extended Data Fig. 3, in which it is clear that the Ramsey contrast decays more quickly as the density of the molecules increases, and, therefore, the dipolar interactions between the molecules is the limiting factor in the coherence time in this system once a spin echo pulse is added.

Each term V(r i - r j ) in equation ( 1) can be calculated numerically. With |↑⟩ and |↓⟩ as the pseudospin-1/2 system, the spin-exchange interac- tion experienced by two molecules with wavefunctions ψ 1 (r i ) and ψ 2 (r j ) on lattice sites r i , r j is given by ∫

where

Here, d = 3.3 D is the permanent body-frame dipole moment 57 , ε 0 is the vacuum permittivity, and ⟨1, -1, 3/2, 3/2 ↓⟩ is the wavefunction overlap between the target |↓⟩ state in the N = 1 manifold and the |1, -1, 3/2, 3/2⟩ state, leading to J/h = 610 Hz at 60 G and J/h = 659 Hz at 4.1 G. The dominant source of uncertainty in J/h is a difference of approximately 0.1 D between theoretical 57 and experimental 42 values of d, leading to corrections of approximately 6% in J/h for the isotropic and anisotropic cases. Uncertainties in a lat lead to corrections of less than 3%. The molecules occupy the ground centre-of-mass state of their respective lattice sites, so their wavefunctions are approximated by the 3D harmonic oscillator ground state wavefunctions given by ψ(r i ) = ψ(x i )ψ(y i )ψ(z i ). The axial (z) and radial (x, y) trap frequencies are 2 kHz and 9 kHz, respectively. The resulting calculations are shown in Extended Data Fig. 4 and in Fig. 4.

Theoretical data in the figures were generated by using exact numerical diagonalization of the XY and XXZ Hamiltonians for samples of N = 12 spins randomly placed on a square lattice of size L × L with periodic boundaries. The spin-spin interactions for a given displacement are set as described in the previous section, such that the large r form of the Hamiltonians match onto equations ( 1) and ( 2) with | J|/h = 600 Hz. To accommodate long-range interactions and periodic boundaries, the interaction between two spins is set according to the shortest displacement between the spins. Each data point shown for comparison to theory is generated from 2,500 samples. These samples are meant to model local patches of the much larger experimental lattice. Owing to the experimental run times being moderate in units of the interaction time for typical spins, simulating a larger number of spins is not necessary. Integer L is chosen for a sample such that N/L 2 is closest to a target density.

The density profile of molecules in the experimental lattice is non-uniform. To model this in the simulations, we let the density of our samples be a random variable. Lower (higher) density samples represent patches of the experimental lattice that are farther from (closer to) the centre. We assume a radial density profile for the experiment of the form

where ρ max is the peak filling at the centre of the lattice and r max is the distance from the centre at which the density becomes negligible. None of our theoretical results actually depend on r max , and the peak densi- ties for each experiment are given in the main text. Sampling patches of N spins with density ρ from such a density profile corresponds to sampling a radial location according to the probability density p(r) ∝ rρ(r). The factor of r is because the amount of lattice at radius r is proportional to r, and the factor of ρ(r) is because a patch with a predetermined number of spins N covers an area proportional to 1/ρ(r). Now, instead of sampling the position of the patch in the lattice, we change the variables and sample the corresponding density of the patch. This yields a probability density p(ρ) ∝ ρ for ρ ρ ≤ max . Correlation functions are computed in a way that is similar to the analysis of experimental data: for the experimental data we first compute the quantum-and-disorder-averaged values of r r a n n ↑ + ↑ and r n ↑ , then form the position-dependent correlation function

and average that over the lattice. The lattice average is performed last so as to try to delay mixing data from different densities until after the connected correlation function is formed. Similarly, for the simulation data we first compute the connected correlation function independently for groups of samples with the same density, and then we average that density-dependent correlation function, with the weight in the average given by the total lattice area of the samples at that density.

Both experimental and simulated correlation functions are scaled in the same way, namely by dividing the amplitude of the correlations by the lattice-averaged value of the squared density ρ ⟨ ⟩ 2 lat

. As discussed in the main text, for an isolated pair of molecules with displacement a the probability of detecting the molecules in state |↑↑⟩ after a Ramsey sequence with evolution time T is V T ħ sin ( ( ) /4 ) 2 a

. When generalizing to larger systems of molecules, owing to the low filling fractions used in the Ramsey experiments we can approximate our system as consisting of single molecules and isolated pairs. In addition, at the end of our Ramsey sequence, single molecules that do not experience many-body interactions are in |↓⟩ , so they do not contribute to the correlation function. . For low fillings the second term can be neglected, so we normalize our correlations by multiplying it by ρ (⟨ ⟩ )

2 lat -1

. For the numerical simulations with a density given by equation (7) , which allows for convergence with fewer experimental repetitions.

A single-fit parameter is introduced to scale the amplitude of the numerically simulated correlations to the experimental data. The amplitude scale factors are obtained by fitting the experimental data for the displacements in the line plots simultaneously and are 0.59(2) for Fig. 3, 0.56(3) for Fig. 4b-d, 0.33(3) for Fig. 4f-h and 0.83 (5) for Fig. 5. The deviation of the experimental correlation amplitudes from the numerical simulations may be due to several factors. The correlation data were acquired over several hours for each set (greater than 24 hours for certain sets). Because the normalization factors were typically determined by imaging the density distribution of the molecules at the beginning of each data set, slow downwards drifts in the molecule number during a data set could affect the amplitude of the correlations. For example, the lattice-averaged correlation before density normalization scales as ρ ⟨ ⟩ 2 lat , so if the density during the correlation data set drifts downward by 30% from the initial density, then the correlation amplitude would drift downward by approximately 50%. An additional source of deviation between the experimental and numerically simulated correlation amplitudes is the possibility of clustering of molecules arising from the detailed spatial overlap of the atomic clouds before molecule formation. This would lead to a modified scaling of the correlation amplitude with density, which would not be captured by the numerical simulations. In the future, the influence of density drifts and molecule clustering could be largely eliminated by implementing spin-resolved imaging. This would enable the post-selection of isolated pairs of molecules in each experimental realization, simplifying the study of the dynamics of the molecular spin system.

In Fig. 5 we report measurements of the dynamics of an effective spin-anisotropic Heisenberg model engineered by using Floquet driving. However, in general there are several ways that the effective Floquet Hamiltonian can deviate from the desired target Hamiltonian over the evolution time of the experiment. Possible errors include the Floquet period not being sufficiently small compared to the fastest coupling timescale in the original Hamiltonian, or if the Bloch sphere rotations are too slow or have improper timing. To check this, we include additional numerical simulations comparing the target XXZ Hamiltonian to different Floquet driving protocols. Extended Data Fig. 5a,b compares the simulated dynamics of the XXZ Hamiltonian to a Floquet drive protocol with the same cycle duration as the experiment from Fig. 5 (Floquet period T F = 1.1124 ms) but a much faster π-pulse duration of 5.8 μs. The simulated dynamics of the magnetization and correlation oscillation show excellent agreement, indicating that over the timescale of the experiment, the Floquet period is sufficiently short to yield a good approximation to the target XXZ Hamiltonian.

Next, in Extended Data Fig. 5c,d we compare the dynamics of the target XXZ Hamiltonian to a Floquet drive with the same cycle length but a more realistic pulse duration of 58 μs. Here, we see slight deviations between the target XXZ Hamiltonian and the Floquet drive in the nearest-neighbour correlation oscillations, which are also present in the experimental data from Fig. 5c. These simulations indicate that the finite pulse duration is more likely to be the cause for this discrepancy than the Floquet cycle being insufficiently fast. Therefore, in future work, it would be valuable to perform the same Floquet experiment with varying microwave Rabi frequencies to test its effect on the Floquet Hamiltonian.

The experimental data from Fig. 5 also shows a slight oscillation in the magnetization of the | X + ⟩ state, which should be a conserved quan- tity in the target Hamiltonian in equation ( 2). Another possible source of error for the experiment that could potentially explain this effect is imperfect π-pulse timing, leading to inadvertent mixing of the spin components and deviations from the target Hamiltonian. We test this for each data set by repeatedly driving transitions between ↑⟩ | and ↓⟩ | and measuring the number of molecules in |↑⟩ after N π pulses. For example, we show in Extended Data Fig. 6 the results of this measurement at 60 G with square microwave pulses of duration 58.3 μs separated by an interval time τ = 100 μs. These results indicate that for the maximum number of 25 π pulses used in the Floquet experiments from Fig. 5, the pulse timing error was negligible, and, therefore, cannot explain any slight deviations from the target Hamiltonian. Future experiments will be required to more fully characterize the precise manner in which the Floquet drive breaks down as an approximation to the target Hamiltonian. For example, by probing the Floquet dynamics to much longer times and varying the duration of the pulses as well as the total Floquet period, it should be possible to observe a more dramatic breakdown of the approximation to the target Hamiltonian to pinpoint the limitations to this technique.