Controlling Dipolar Exchange Interactions in a Dense 3D Array of Large Spin Fermions

TL;DR

This study demonstrates control over dipolar exchange interactions in a dense 3D array of large-spin fermionic erbium atoms, enabling tunable quantum spin dynamics and potential applications in quantum information and exotic magnetic phases.

Contribution

We experimentally control dipolar spin-exchange interactions in a dense 3D lattice of erbium atoms by tuning dipole orientation and initial states, with results matching advanced numerical models.

Findings

Controlled spin dynamics via dipole orientation and initial state tuning

Ability to switch dipolar interactions on and off rapidly

Excellent agreement between experiment and numerical simulations

Abstract

Dipolar interactions are ubiquitous in nature and rule the behavior of a broad range of systems spanning from energy transfer in biological systems to quantum magnetism. Here, we study magnetization-conserving dipolar induced spin-exchange dynamics in dense arrays of fermionic erbium atoms confined in a deep three-dimensional lattice. Harnessing the special atomic properties of erbium, we demonstrate control over the spin dynamics by tuning the dipole orientation and changing the initial spin state within the large 20 spin hyperfine manifold. Furthermore, we demonstrate the capability to quickly turn on and off the dipolar exchange dynamics via optical control. The experimental observations are in excellent quantitative agreement with numerical calculations based on discrete phase-space methods, which capture entanglement and beyond-mean field effects. Our experiment sets the stage for future explorations of rich magnetic behaviors in long-range interacting dipoles, including exotic phases of matter and applications for quantum information processing.

Tables

Table 1. Table 1 : Percentage of lost atoms and extracted, effective filling fraction ν 𝜈 \nu for different | m F ⟩ ket subscript 𝑚 𝐹 |m_{F}\rangle states at t = t fit 𝑡 subscript 𝑡 fit t=t_{\rm{fit}} and t = 100 𝑡 100 t=100\, ms.

$m_F^0$	$t_{\mathrm{fit}}$ (ms)	$N_{\mathrm{loss}}\left(t_{\mathrm{fit}}\right)$ (%)	$\nu$ $\left(0\right)$	$\nu$ $\left(t_{\mathrm{fit}}\right)$	$N_{\mathrm{loss}}(100$ms) (%)	$\nu$ $(100$ms)
$-\frac{17}{2}$	$34.2$	$1.8$	$0.7$	$0.68$	$5.3$	$0.66$
$-\frac{13}{2}$	$15.7$	$7.2$	$0.69$	$0.64$	$19.6$	$0.55$
$-\frac{9}{2}$	$11.3$	$8.7$	$0.67$	$0.62$	$25.1$	$0.51$
$-\frac{5}{2}$	$9.6$	$13.7$	$0.66$	$0.58$	$27.7$	$0.48$
$-\frac{1}{2}$	$9.0$	$12.2$	$0.65$	$0.57$	$36.1$	$0.42$
$\frac{1}{2}$	$9.0$	$13.5$	$0.65$	$0.56$	$36.7$	$0.41$
$\frac{3}{2}$	$9.2$	$9.2$	$0.64$	$0.58$	$34.0$	$0.43$
$\frac{9}{2}$	$11.3$	$6.7$	$0.63$	$0.59$	$23.1$	$0.49$
$\frac{13}{2}$	$15.7$	$4.9$	$0.62$	$0.59$	$21.9$	$0.48$

Table 2. Table 2 : Percentage of lost atoms and extracted effective filling fraction ν 𝜈 \nu for different Θ Θ \Theta at t = t fit 𝑡 subscript 𝑡 fit t=t_{\rm{fit}} and t = 100 𝑡 100 t=100\, ms.

$\Theta$ (∘)	$t_{\mathrm{fit}}$ (ms)	$N_{\mathrm{loss}}\left(t_{\mathrm{fit}}\right)$ (%)	$\nu$ $\left(0\right)$	$\nu$ $\left(t_{\mathrm{fit}}\right)$	$N_{\mathrm{loss}}(100$ms) (%)	$\nu$ $(100$ms)
$0$	$26.8$	$13.8$	$0.7$	$0.60$	$30.1$	$0.49$
$10$	$30.1$	$11.9$	$0.7$	$0.61$	$25.6$	$0.52$
$20$	$36.7$	$6.9$	$0.7$	$0.65$	$18.2$	$0.57$
$30$	$47.0$	$8.6$	$0.7$	$0.64$	$17.8$	$0.57$
$35$	$52.0$	$6.7$	$0.7$	$0.65$	$11.4$	$0.62$
$40$	$53.2$	$7.2$	$0.7$	$0.65$	$13.2$	$0.60$
$50$	$46.1$	$12$	$0.7$	$0.61$	$17.2$	$0.58$
$60$	$37.0$	$8.9$	$0.7$	$0.63$	$19.2$	$0.56$
$80$	$30.0$	$10.4$	$0.7$	$0.62$	$19.3$	$0.56$

Equations

$$\hat{H}$$

$$n_{m_F^0}(t)$$

$$U_t=\frac{3m_F^2-F(F+1)}{F(2F-1)}\frac{3{\cos }^2{\theta }_p-1}{2}{\alpha }_t(\omega ),$$

$$\hat{H}$$

$$V_{\mathbf{i},\mathbf{j}}$$

$$V_{dd}$$

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

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

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

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

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

$$c_{\mu }^i$$

$$i\hslash \frac{d{\Lambda }_{\mu }^i}{dt}$$

$${\lambda }_{\mu }^i(t)$$

$$n_{m_F}(t)$$

$$n_{m_F^0}(t)$$

$$V_{\mathrm{eff}}^2$$

$$\gamma (m_F^0)$$

$$n_{m_F^0}^{\mathrm{Mean}-\mathrm{Field}}(t)$$

Full text

Spin lattice models of localized magnetic moments (spins), which interact with one another via exchange interactions, are paradigmatic examples of strongly correlated many-body quantum systems. Their implementation in clean, isolated, and fully controllable lattice confined ultra-cold atoms opens a path for a new generation of synthetic quantum magnets, featuring highly entangled states, especially when driven out-of-equilibrium, with broad applications ranging from precision sensing and navigation, to quantum simulation and quantum information processing Bloch (2008); Gross and Bloch (2017). However, the extremely small energy scales associated with the nearest-neighbor spin interactions in lattice-confined atoms with dominant contact interactions, have made the observation of quantum magnetic behaviors extremely challenging Bloch et al. (2008); Greif et al. (2013). On the contrary, even under frozen motional conditions, dipolar gases, featuring long-range and anisotropic interactions, offer the opportunity to bring ultra-cold systems several steps ahead towards the ambitious attempt of modeling and understanding quantum magnetism. While great progress in studying quantum magnetism has been achieved using arrays of Rydberg atoms Zeiher et al. (2017); Bernien et al. (2017); Barredo et al. (2018); Guardado-Sanchez et al. (2018), trapped ions Neyenhuis et al. (2017); Blatt and Roos (2012); Britton et al. (2012), polar molecules Yan et al. (2013); Hazzard et al. (2014), and spin-$3$ bosonic chromium atoms de Paz et al. (2013, 2016); Lepoutre et al. (2018), most of the studies so far have been limited to spin-$1/2$ mesoscopic arrays of at the most few hundred particles or to macroscopic but dilute ($<0.1$ filling fractions) samples of molecules in lattices.

In this work, we report the first investigations of non-equilibrium quantum magnetism in a dense array of fermionic magnetic atoms confined in a deep three-dimensional optical lattice. Our platform realizes a quantum simulator of the long-range XXZ Heisenberg model. The simulator roots on the special atomic properties of 167Er, whose ground state bears large angular momentum quantum numbers with $I=7/2$ for the nuclear spin and $J=6$ for the electronic angular momentum, resulting in a $F=19/2$ hyperfine manifold, as depicted in Fig. 1A. Such a complexity enables new control knobs for quantum simulations. First, it is responsible for the large magnetic moment in Er. Second, it gives access to a fully controllable landscape of $20$ internal levels, all coupled by strong magnetic dipolar interactions, up to $49$ times larger than the ones felt by $F=1/2$ alkali atoms in the same lattice potential Stamper-Kurn and Ueda (2013). Finally, it allows fast optical control of the energy splitting between the internal states which can be tuned on and off resonance using the large tensorial light shift Becher et al. (2018), which adds to the usual quadratic Zeeman shift.

Using all these control knobs, we explore the dipolar exchange dynamics and benchmark our simulator with an advanced theoretical model, which takes entanglement and beyond mean-field effects into account sup . In particular, we initialize the system into a desired spin state and activate the spin dynamics, while the motional degree of freedom mainly remains frozen. Here, we study the spreading of the spin population in the different magnetic sublevels as a function of both the specific initial spin state and the dipole orientation. We demonstrate that the spin dynamics at short evolution time follows a scaling that is invariant under internal state initialization (choice of macroscopically populated initial Zeeman level) and that is set by the effective strength of the dipolar coupling. On the contrary, at longer times, the many-body dynamics is affected by the accessible spin space and the long-range character of dipolar interactions beyond nearest neighbors. Finally, we show temporal control of the exchange dynamics using off resonant laser light.

The XXZ Heisenberg model that rules the magnetization-conserving spin dynamics of our system can be conveniently written using spin-$19/2$ dimensionless angular momentum operators $\hat{\mathbf{F}}_{i}=\{\hat{F}^{x}_{i},\hat{F}^{y}_{i},\hat{F}^{z}_{i}\}$, acting on site $i$ and satisfying the commutation relation $[\hat{F}_{i}^{x},\hat{F}_{i}^{y}]=i\hat{F}_{i}^{z}$. We use the eigenbasis of $\hat{F}^{z}$ denoted as $|m_{F}\rangle$ with $0\leq|m_{F}|\leq F$  Auerbach (1994); Dutta et al. (2015); sup :

[TABLE]

The coupling constants $V_{i,j}\,=\,V_{dd}d_{y}^{3}\frac{1-3\cos^{2}(\theta_{i,j})}{r_{ij}^{3}}$, describe the direct dipole-dipole interactions (DDI), which have long-range character and thus couple beyond nearest neighbors. The dipolar coupling strength between two dipoles located at $\vec{r}_{i}$ and $\vec{r}_{j}$ depends on their relative distance $r_{ij}=|\vec{r}_{i}-\vec{r}_{j}|$ and on their orientation, described by the angle $\theta_{i,j}$ between the dipolar axis, set by the external magnetic field, and the interparticle axis; see Fig. 1B. Here, $V_{dd}\,\approx\,\frac{\mu_{0}g_{F}^{2}\mu_{B}^{2}}{4\pi d_{y}^{3}}$ denotes the dipolar coupling strength, with $g_{F}\approx 0.735$ for 167Er, $\mu_{0}$ the magnetic permeability of vacuum, $\mu_{B}$ the Bohr magneton, and $d_{y}$ the shortest lattice constant. The $\hat{F}_{i}^{z}\hat{F}_{j}^{z}$ terms in the Hamiltonian account for the diagonal part of the interactions while the $\hat{F}^{+}_{i}\hat{F}^{-}_{j}+\hat{F}^{-}_{i}\hat{F}^{+}_{j}$ terms describe dipolar exchange processes. The second sum denotes the single particle quadratic term $\delta_{i}(\hat{F}^{z}_{i})^{2}$ with $\delta_{i}=\delta_{i}^{Z}+\delta_{i}^{T}$, accounting for the quadratic Zeeman effect $\propto\delta_{i}^{Z}$ and tensorial light shifts $\propto\delta_{i}^{T}$. These two contributions can be independently controlled in our experiment.

The quadratic Zeeman shift allows us to selectively prepare all atoms in one target state of the spin manifold sup . The tensorial light shift can compete or cooperate with the quadratic Zeeman shift and can be used as an additional control knob to activate/deactivate the exchange processes. Note that, for all measurements, a large linear Zeeman shift is always present, but since it does not influence the spin-conserving dynamics, it is omitted in Eq. 1.

In the experiment, we first load a spin-polarized quantum degenerate Fermi gas of $\approx 10^{4}$ Er atoms into a deep 3D optical lattice, following the scheme of Ref. Baier et al. (2018). The cuboid lattice geometry with lattice constants $(d_{x},d_{y},d_{z})=(272,266,544)\,$nm results in weakly coupled 2D planes, with typical tunneling rates of $\sim 10\,$Hz inside the planes and $\sim\,$mHz between them sup . The external magnetic field orientation, setting the quantization axis as well as the dipolar coupling strengths, is defined by the polar angles $\Theta$ and $\phi$ in the laboratory frame; see Fig. 1B. The fermionic statistics of the atoms enables us to prepare a dense band insulator with at most one atom per lattice site, as required for a clean implementation of the XXZ Heisenberg model. This is an advantage of fermionic atoms as compared to bosonic systems, which typically require filtering protocols to remove doublons Lepoutre et al. (2018).

Our experimental sequence to study the spin dynamics is illustrated in Fig. 1C. In particular, we prepare the system into the targeted $m_{F}^{0}$ state by using the lattice-protection protocol demonstrated in Ref. Baier et al. (2018); see also Ref. sup . At the end of the preparation, the majority of atoms are in the desired $m_{F}^{0}$ ($>80\%$) at $B\approx 4\,$G. We note that atom losses during the spin preparation stage reduces the filling factor to about 60% of the initial one sup . We then activate the spin dynamics by quenching the magnetic field to a value for which $\bar{\delta}=\sum_{i}\delta_{i}=0$, providing a resonance condition for the magnetization-conserving spin-exchange processes; see Fig. 2A. After a desired time of evolution, we stop the dynamics by rapidly increasing the magnetic field, leaving the resonance condition. We finally extract the atom number in each spin state via a spin-resolved band-mapping technique sup and derive the relative state populations by normalization to the initial total atom number.

We now probe the evolution of the spin-state population as a function of the hold time on resonance. We observe a redistribution of the population from the initial state to multiple neighboring states in $m_{F}$ space, as exemplary shown for an initial state of $|\text{--}13/2\rangle$ in Fig. 2B-C. The dynamics preserves the total magnetization; see inset of Fig. 2B. We observe similar behavior independently of the initialized $m_{F}^{0}$ states. The spin transfer happens sequentially. At short times it is dominated by the transfer to states directly coupled by the dipolar exchange Hamiltonian, i. e. those ones which differ by plus/minus one unit of angular momentum ($\Delta m_{F}=\pm 1$). At longer times, subsequent processes transfer atoms to states with $|\Delta m_{F}|\geq 2$; see Fig. 2C-D.

To benchmark our quantum simulator, we use a semiclassical phase-space sampling method, the so called generalized discrete truncated Wigner approximation (GDTWA) Lepoutre et al. (2018); Schachenmayer et al. (2015a, b); Polkovnikov (2010); sup . The method accounts for quantum correlation in the many-body dynamics and is adapted to tackle the complex dynamics of a large-spin system in a regime where exact diagonalization techniques are impossible to implement with current computers. The GDTWA calculations take into account actual experimental parameters such as spatial inhomogeneites, density distribution after the lattice loading, initial spin distribution, and effective lattice filling, including the loss during the spin preparation protocol sup . Figure 2B shows the experimental dynamics together with the GDTWA simulations. Although the model does not include corrections due to losses and tunneling during the dynamics, it successfully captures the behavior of our dense system not only at short time, but also at long time, where the population dynamics slows down and starts to reach an equilibrium. Similar level of agreement between experiment and theory is shown in Fig. 2C-D where we directly compare the spreading of the spin population as a function of time.

The important role of quantum effects in the observed spin dynamics can be illustrated by contrasting the GDTWA simulation with a mean-field calculation. Indeed, the mean-field calculation fails in capturing the system behavior. It predicts a too slow population dynamics for non-perfect spin-state initialization, as in the experiments shown in Fig. 2B, and no dynamics for the ideal case where all atoms are prepared in the same internal state sup . To emphasize the beyond nearest-neighbor effects, we also compare the experiment with a numerical simulation that only includes nearest-neighbor interactions (NNI-GDTWA). Here, we again observe a very slow spin evolution, which largely deviates from the measurements. The agreement of the full GDTWA predictions with our experimental observations points to the long-range many-body nature of the underlying time evolution. Our theory calculations also support the built up of entanglement during the observed time evolution.

Different spin configurations feature distinct effective interaction strengths, which also depend on the orientation of the dipoles with respect to the lattice. We demonstrate our ability to control this interaction, which governs the rate of population exchange, by the choice of the initial $m_{F}^{0}$ state and the orientation of the external magnetic field. This capability provides us with two tuning knobs to manipulate dipolar exchange interactions in our quantum simulator. Figure 3A-F plots the dynamics of the populations for three neighboring spin states after the quench, starting from different initial spin states. Solid lines show the results of the full-GDTWA calculations. For each initial $m_{F}^{0}$, we find a remarkable agreement between theory and experiment. We observe a strong speedup for states with large spin projections perpendicular to the quantization axis, as it is expected from the expectation value of $\hat{F}^{+}_{i}\hat{F}^{-}_{j}$, which gives a prefactor $\gamma(m_{F}^{0})=\sqrt{F(F+1)-m_{F}^{0}(m_{F}^{0}+1)}\sqrt{F(F+1)-m_{F}^{0}(m_{F}^{0}-1)}$. The initial dynamics can be well described by a perturbative expansion up to the second order sup , resulting in the analytic expression for the normalized population $n_{m_{F}}(t)$ of the initial state:

[TABLE]

Here, $V^{2}_{\text{eff}}\equiv\frac{\gamma^{2}(m_{F}^{0})}{8N}\sum_{i,j\neq i}V_{ij}^{2}$ is the overall effective interaction strength summed over $N$ atoms and $n_{m_{F}^{0}}(0)$ denotes the purity of the initial state preparation. For a quantitative analysis of the early-time spin evolution, we compare the theoretically calculated $V_{\text{eff}}$ from the initial atomic distribution used in the GDTWA model with the one extracted from a fit of Eq. 2 to the experimental data. Here we consider the data up to $t<0.5\frac{\hbar}{V_{\text{eff}}}$ estimated using the theoretically calculated $V_{\text{eff}}$ Not . Figure 3G plots both, the theoretical and experimental $V_{\text{eff}}$ as a function of the initial $m_{F}^{0}$ and highlights once more their quantitative agreement. The interaction parameter $V_{\text{eff}}$, can also be used to rescale the time axis. As shown in Fig. 3H, all data sets now collapse onto each other for $\frac{tV_{\text{eff}}}{\hbar}<0.5$, revealing the invariant character of the short-time dynamics under the internal state initialization. At longer timescales, the theory shows that the timetraces start to deviate from each others and saturate to different values, indicating that thermal-like states are on reach. In the experiment, we observe a similar behavior but here the saturation value might also be affected by losses and residual tunneling.

Because of the anisotropic character of the DDI, the strength of the dipolar exchange can be controlled by changing the angle $\Theta$; see Fig. 1B. As exemplary shown in Fig. 4A for $|\text{--}17/2\rangle$, the observed evolution speed of the spin populations strongly depends on $\Theta$, changing by about a factor of $2$ between $\Theta=40^{\circ}$ and $80^{\circ}$. The GDTWA results show a very good quantitative agreement with the experiment. We repeat the above measurements for different values of $\Theta$ and we extract $V_{\text{eff}}$; Fig. 4B. It is worth to notice that, while the dipolar interactions can be completely switched off at a given angle in a 1D chain, in a 3D system the situation is more complicated. However, as expected by geometrical arguments, we observe that the overall exchange strength becomes minimal for a specific dipole orientation $(\Theta_{c}\approx 35^{\circ},\phi_{c}=45^{\circ})$. We compare our measured $V_{\text{eff}}$ with the ones calculated from the initial spin distribution, which is a good quantity to describe the early time dynamics. Theory and experiment show a similar trend, in particular reaching a minimum at about $\Theta_{c}$. Note that the simple analytic formula (Eq. 2), used for fitting the data, deviates from the actual evolution at longer times. This leads to a small down-shift of the experimental values sup .

Finally, we demonstrate fast optical control of the spin dynamics relying on the remarkably large tensorial light shift of erbium compared to alkali atoms. As shown in Fig. 4C, we can almost fully suppress the spin exchange dynamics by suddenly switching on a homogeneous light field after an initial evolution time on resonance. Therefore the tensorial light shift, inducing a detuning from the resonance condition (see inset), allows a full spatial and temporal control over the exchange processes as the light power can be changed orders of magnitude faster than the typical interaction times and can address even single lattice sites in quantum gas microscopes. This capability can be an excellent resource for quantum information processing, i. e. we could use dipolar exchange processes to efficiently prepare highly entangled states between different parts of a quantum system and then store the quantum information at longer times by turning the interactions off.

The excellent agreement between the experiment and the theory, not only verifies our quantum simulator but sets the stage towards its use for the study of new regimes intractable to theory. For example by operating at weaker lattice when motion is involved, the dynamics is no longer described by a spin model but by a high spin Fermi-Hubbard model with long-range interactions. Alternatively by treating the internal hyperfine levels as a synthetic dimension Mancini et al. (2015) while adding Raman transitions to couple them, one could engineer non-trivial synthetic gauge field models even when tunneling is only allowed in one direction. Moreover, the demonstrated control over the different hyperfine level structure, our capability to tune the strength of the direct dipolar exchange coupling via the magnetic field angle, and the possibility of the dynamical and spatial control of the resonance condition via tensorial light shifts make our system a potential resource for quantum information processing with a qudit-type multi-level encoding using the 20 different interconnected hyperfine levels Lanyon et al. (2008); Brion et al. (2007); Kues et al. (2017).

• Correspondence and requests for materials should be addressed to M. J. M. (email: manfred.mark@uibk.ac.at).

References and Notes

• Bloch (2008) I. Bloch, Nature 453, 1016 (2008).
• Gross and Bloch (2017) C. Gross and I. Bloch, Science 357, 995 (2017).
• Bloch et al. (2008) I. Bloch, J. Dalibard,  and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
• Greif et al. (2013) D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell,  and T. Esslinger, Science 340, 1307 (2013).
• Zeiher et al. (2017) J. Zeiher, J.-y. Choi, A. Rubio-Abadal, T. Pohl, R. van Bijnen, I. Bloch,  and C. Gross, Phys. Rev. X 7, 041063 (2017).
• Bernien et al. (2017) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić,  and M. D. Lukin, Nature 551, 579 (2017).
• Barredo et al. (2018) D. Barredo, V. Lienhard, S. de Leseleuc, T. Lahaye,  and A. Browaeys, Nature 561, 79 (2018).
• Guardado-Sanchez et al. (2018) E. Guardado-Sanchez, P. T. Brown, D. Mitra, T. Devakul, D. A. Huse, P. Schauß,  and W. S. Bakr, Phys. Rev. X 8, 021069 (2018).
• Neyenhuis et al. (2017) B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C. Lee, P. Richerme, Z.-X. Gong, A. V. Gorshkov,  and C. Monroe, Science Advances 3, e1700672 (2017).
• Blatt and Roos (2012) R. Blatt and C. F. Roos, Nature Physics 8, 277 (2012).
• Britton et al. (2012) J. W. Britton, B. C. Sawyer, A. C. Keith, C. C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk,  and J. J. Bollinger, Nature 484, 489 (2012).
• Yan et al. (2013) B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin,  and J. Ye, Nature 501, 521 (2013).
• Hazzard et al. (2014) K. R. A. Hazzard, B. Gadway, M. Foss-Feig, B. Yan, S. A. Moses, J. P. Covey, N. Y. Yao, M. D. Lukin, J. Ye, D. S. Jin,  and A. M. Rey, Phys. Rev. Lett. 113, 195302 (2014).
• de Paz et al. (2013) A. de Paz, A. Sharma, A. Chotia, E. Maréchal, J. H. Huckans, P. Pedri, L. Santos, O. Gorceix, L. Vernac,  and B. Laburthe-Tolra, Phys. Rev. Lett. 111, 185305 (2013).
• de Paz et al. (2016) A. de Paz, P. Pedri, A. Sharma, M. Efremov, B. Naylor, O. Gorceix, E. Maréchal, L. Vernac,  and B. Laburthe-Tolra, Phys. Rev. A 93, 021603 (2016).
• Lepoutre et al. (2018) S. Lepoutre, J. Schachenmayer, L. Gabardos, B. Zhu, B. Naylor, E. Marechal, O. Gorceix, A. M. Rey, L. Vernac,  and B. Laburthe-Tolra, arxiv 1803.02628 (2018).
• Stamper-Kurn and Ueda (2013) D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. 85, 1191 (2013).
• Becher et al. (2018) J. H. Becher, S. Baier, K. Aikawa, M. Lepers, J.-F. Wyart, O. Dulieu,  and F. Ferlaino, Phys. Rev. A 97, 012509 (2018).
• (19) Materials, methods, and additional theoretical background are available as supplementary material on Science Online.
• Auerbach (1994) A. Auerbach, Interacting Electrons and Quantum Magnetism (Springer-Verlag, New York, 1994).
• Dutta et al. (2015) O. Dutta, M. Gajda, P. Hauke, M. Lewenstein, D.-S. Lühmann, B. A. Malomed, T. Sowiński,  and J. Zakrzewski, Reports on Progress in Physics 78, 066001 (2015).
• Baier et al. (2018) S. Baier, D. Petter, J. H. Becher, A. Patscheider, G. Natale, L. Chomaz, M. J. Mark,  and F. Ferlaino, Phys. Rev. Lett. 121, 093602 (2018).
• Schachenmayer et al. (2015a) J. Schachenmayer, A. Pikovski,  and A. M. Rey, New Journal of Physics 17, 065009 (2015a).
• Schachenmayer et al. (2015b) J. Schachenmayer, A. Pikovski,  and A. M. Rey, Phys. Rev. X 5, 011022 (2015b).
• Polkovnikov (2010) A. Polkovnikov, Annals of Physics 325, 1790 (2010).
• (26) We note that dissipationless dynamics is the most justified on this time scale as losses stay below 14%.
• Mancini et al. (2015) M. Mancini, G. Pagano, G. Cappellini, L. Livi, M. Rider, J. Catani, C. Sias, P. Zoller, M. Inguscio, M. Dalmonte,  and L. Fallani, Science 349, 1510 (2015).
• Lanyon et al. (2008) B. P. Lanyon, M. Barbieri, M. P. Almeida, T. Jennewein, T. C. Ralph, K. J. Resch, G. J. Pryde, J. L. O’Brien, A. Gilchrist,  and A. G. White, Nature Physics 5, 134 (2008).
• Brion et al. (2007) E. Brion, K. Mølmer,  and M. Saffman, Phys. Rev. Lett. 99, 260501 (2007).
• Kues et al. (2017) M. Kues, C. Reimer, P. Roztocki, L. R. Cortés, S. Sciara, B. Wetzel, Y. Zhang, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, J. Azaña,  and R. Morandotti, Nature 546, 622 (2017).
• Aikawa et al. (2014) K. Aikawa, A. Frisch, M. Mark, S. Baier, R. Grimm,  and F. Ferlaino, Phys. Rev. Lett. 112, 010404 (2014).
• Baier et al. (2016) S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Z. Cai, M. Baranov, P. Zoller,  and F. Ferlaino, Science 352, 201 (2016).
• Bertlmann and Krammer (2008) R. A. Bertlmann and P. Krammer, Journal of Physics A: Mathematical and Theoretical 41, 235303 (2008).

Supplementary Materials

S1 Experimental setup and lattice loading

In our experiment we start with a degenerate Fermi gas of about $2.4\times 10^{4}$ 167Er atoms in the lowest spin state $|F=19/2,m_{F}=\text{--}19/2\rangle=|\text{--}19/2\rangle$ and a temperature of $T\approx 0.3\,T_{\text{F}}$ Aikawa et al. (2014); Baier et al. (2018). The atoms are confined in a crossed optical dipole trap (ODT) and the trap frequencies are $(\nu_{\perp},\nu_{\parallel},\nu_{z})=(63(1),36(2),137(1))\,\rm Hz$, where $\nu_{\perp}$ ($\nu_{\parallel}$) are the trap frequencies perpendicular to (along) the horizontal ODT and $\nu_{z}$ is measured along the vertical direction defined by gravity. We load the atomic sample adiabatically into a 3D lattice that is created by two retro-reflected laser beams at $532\,$nm in the x-y plane and one retro-reflected laser beam at $1064\,$nm nearly along the z direction, defined by gravity and orthogonal to the x-y plane. Note, that due to a small angle of about $10^{\circ}$ between the vertical lattice beam and the z axis we obtain a 3D-lattice, slightly deviating from an ideal situation of a rectangular unit cell and our parallelepipedic cell has the unit lattice distances of $d_{x}\,=\,272\,{\rm nm}$, $d_{y}\,=\,266\,{\rm nm}$, and $d_{z}\,=\,544\,{\rm nm}$. The lattice geometry is similar to the one used in our previous works Baier et al. (2016, 2018). We ramp up the lattice beams in $150\,$ms to their final power and switch off the ODT subsequently in $10\,$ms and wait for $500\,$ms to eliminate residual atoms in higher bands Baier et al. (2018). For our typical lattice depths used in the measurements reported here of $(s_{x},s_{y},s_{z})\,=\,(20,20,80)$, where $s_{i}$ with $i\in{x,y,z}$ is given in the respective recoil energy $\mathrm{E}_{\mathrm{R},i}$ with $E_{\text{R};x,y}=h\times 4.2\,\rm kHz$ and $E_{\text{R};z}=h\times 1.05\,\rm kHz$, the atoms can be considered pinned on single lattice sites with low tunneling rates $J_{x,y}=h\times 10.5\,\rm Hz$ and $J_{z}=h\times 1\,\rm mHz$.

S2 State preparation and detection efficiency

To prepare the atoms in the desired spin state, after loading them into the lattice, we use a technique based on a rapid-adiabatic passage (RAP). During the full preparation procedure, the optical lattice operates as a protection shield to avoid atom loss and heating due to the large density of Feshbach resonances Baier et al. (2018). To reach a reliable preparation with high fidelity it is necessary that the change in the energy difference between subsequent neighboring spin states is sufficiently large. Therefore, we ramp the magnetic field in $40\,\rm ms$ to $40.6\,\rm G$ to get a large enough differential quadratic Zeeman shift, which is on the order of $\approx h\times 40\,$kHz. After the magnetic field ramp we wait for $80\,\rm ms$ to allow the latter to stabilize before performing the RAP procedure. The follow up RAP protocol depends on the target state. For example, to transfer the atoms from $|\text{--}19/2\rangle$ into the $|\text{--}7/2\rangle$ state, we apply a radio-frequency (RF) pulse at $41.31\,\rm MHz$ and reduce the magnetic field with a linear ramp, e. g. by $500\,\rm mG$ in $42\,\rm ms$. The variation of the magnetic field is analogous to the more conventional scheme where the frequency of the RF is varied (see Fig. S1A). For the preparation of higher (lower) spin states we perform a larger (smaller) reduction of the magnetic field on a longer (shorter) timescale. After the RAP ramp we switch off the RF field and ramp the magnetic field in $10\,\rm ms$ to $B=3.99\,\rm G$. Here we wait again for $100\,$ms to allow the magnetic field to stabilize. During the ramp up and down to $40\,$G of the magnetic field we loose $25(2)\,\%$ of the atoms. We attribute this loss mainly to the dense Feshbach spectra that we are crossing when ramping the magnetic field. The exact loss mechanism has not been yet identified, constituting a topic of interest for latter investigation. At $3.99\,$G, before switching on the spin dynamics, about $1.7\times 10^{4}$ atoms remain in the lattice. The losses affect the lattice filling at which the spin dynamics occur. Our simulations account for this initially reduced filling; see Sec. S8.

Additionally to the losses due to the magnetic field ramps, we also observe losses caused by the RAP itself. To quantify the preparation efficiency, i. e . the loss of atoms due to the spin preparation via RAP as a function of the target $m_{F}$ state, we perform additional measurements where we either do not perform a RAP or where we add an inverse RAP to transfer all atoms back into the $|\text{--}19/2\rangle$ state. By comparing the atom numbers from measurements without RAP and with double-RAP and assuming that the loss process is symmetric, we derive the preparation efficiency as plotted in Fig. S1B. We also account for the difference in the counting efficiency of the individual spin states, which arises from different effective scattering cross sections for the imaging light. Here we compare the measured atom number in a target $m_{F}$ state to the expected atom number taking into account the previously determined preparation efficiency as discussed above and the atom number without RAP; see Fig. S1C.

The counting and preparation efficiencies are determined for the $|\text{--}17/2\rangle$, $|\text{--}15/2\rangle$, $|\text{--}9/2\rangle$, and $|9/2\rangle$ states and interpolated assuming a linear dependency of these efficiencies on $m_{F}$ (see Fig. S1B,C). We estimate the preparation efficiency of the respective mF state to be $0.86(1)-0.008(1)~{}\times~{}\text{m}_{F}$. We attribute the lower preparation efficiency for higher spin states to the larger number of avoided crossing between spin states that come into play during the RAP procedure. Overall we expect that the lattice filling over the whole sample, taking into account the losses due to magnetic field ramping and spin state preparation, reduces from close to unity down to a value between $0.6$ and $0.7$; see also Tab. S1 – S2.

S3 Quench protocol and detection sequence

In our experiment we exploit both, the light and the magnetic shifts of the energies of each spin state to reach a resonant condition where the energy difference between neighboring spin states is cancelled and therefore spin changing dynamics preserving the total magnetization become energetically allowed. In particular we exploit the tensorial light shift of the spin states energies Becher et al. (2018)

[TABLE]

present in atomic erbium to initialize the dynamic evolution of the spin population. The tensorial light shift depends quadratically on the $m_{F}$ state as well as on the angle $\theta_{p}$ between the magnetic field axis and the axis of polarization of the light. Here, $\alpha_{t}$ refers to the tensorial polarizability coefficient and $\omega$ to the light frequency. After the preparation of the respective spin state we start all our measurements at $B=3.99\,\rm G$, pointing in the $z$ direction. However, to reach the resonance condition we use two slightly different quench protocols for the measurement sets with fixed $\Theta=0^{\circ}$ for the different initial spin state and for the sets of measurements where $|m_{F}\rangle=|\text{--}17/2\rangle$ and $\Theta\in(0^{\circ},80^{\circ})$. The measurement sequences differ on the one hand by the way we jump on resonance to initialize the spin dynamics and on the other hand by shining in an additional laser beam of wavelength $1064\,\rm nm$ and power of $7\,\rm W$. This additional light is necessary because changing $\Theta$ reduces simultaneously $\theta_{p}$ resulting in a smaller tensorial light shift and therefore in a shift of the resonance position to lower magnetic field values. For large $\Theta$ the light shift of the lattice beams is smaller and therefore the resonance is very close to $0\,$G which we want to avoid to prevent spin relaxation processes. For the sets of measurements where we keep $\Theta=0^{\circ}$ but vary the initial $m_{F}^{0}$ state we quench the magnetic field directly after the preparation, from $3.99\,\rm G$ to resonance. In contrast we use a different approach for the measurements where $\Theta$ is varied. After the preparation we ramp in $10\,\rm ms$ the additional laser beam to $7\,\rm W$. Due to the reduced speed of our magnetic field coils in $x$ and $y$ direction we first rotate the magnetic field such that the transverse components $B_{\rm{x}}$ and $B_{\rm{y}}$ are already at their target values while keeping an additional offset of $2\,$G in the $z$ direction. The quench to resonance is then done using only the coils for the magnetic field in the $z$ direction. The additional offset field of $2\,$G is large enough to suppress dynamics. We measure the evolution of the magnetic field by performing RF spectroscopy and find that for both quench procedures the magnetic field evolves exponentially towards its quench value with $1/e$ decay times of $1.4\,\rm ms$ and $1.2\,\rm ms$, respectively. After holding on resonance for a certain time we quench the magnetic field back to $3.99\,\rm G$ and we rotate the latter back to $\Theta=0^{\circ}$. After a waiting time of $50\,\rm ms$ we perform a band-mapping measurement combined with a Stern-Gerlach technique, i.e. we ramp the lattice down in $1\,\rm ms$ and apply a magnetic field gradient that is large enough to separate the individual spin states after a time of flight (TOF) of $15\,\rm ms$. This allows us to image the first Brillouin-zone for the different spin states. During TOF the magnetic field is rotated towards the imaging axis. We typically record the population of the initially prepared $|m_{F}\rangle$, of its four neighbors, and of $|\text{--}19/2\rangle$ by summing the 2D atomic density over a region of interest. Figure S2 shows examples of the imaging of different spin states for the cases of a non adjusted RAP as well as for the preparation of the atoms in $|\text{--}9/2\rangle$, $|3/2\rangle$, and $|5/2\rangle$. In the case of $|3/2\rangle$ residual atoms in $|\text{--}19/2\rangle$, $|\text{--}17/2\rangle$, and $|5/2\rangle$ are visible due to a non perfect preparation.

S4 Lifetime and losses in the lattice

Off-resonance, i. e.  at a magnetic field of $B=3.99\,$G, we measure the lifetime of the prepared spin state to be on the order of a few seconds, being slightly shorter for higher spin states. Note that, here, we do not observe any population growing in the neighboring spin states. Differently, for the measurements on resonance we observe a faster loss happening on the timescale of the first $20$–$30\,$ms followed by loss at lower speed for the remaining atoms. We fit an exponential decay to extract the atoms loss and change in filling over the timescale that we use to extract $V_{\rm{eff}}$, $t_{\rm{fit}}$, (see S8) as well as over the full $100\,$ms of the dynamics reported in the main text (Fig. 2-3). Table S1 gives the corresponding numbers for the sets of data for the different initial $m_{F}^{0}$ states. During the fitting timescale we observe atoms loss on the order of $5$–$10$%. This atom loss can be converted into a change of the effective filling of the lattice compared to the state obtained after the lattice loading giving a minimum filling of $\nu=0.58$ for the $|m_{F}^{0}\rangle=|1/2\rangle$ case. For longer timescales larger losses in the range between $10-35$% are observed. In general, we note that the amount of loss depends on the initial $m_{F}^{0}$ state, resulting larger for the central $|m_{F}\rangle^{0}$ states. Similar numbers are obtained for the sets of measurement where we vary $\Theta$ (see Tab. S2). The exact mechanism leading to these losses is not yet understood and will be the topic of future studies. Thanks to their limited importance over the early time dynamics, we here compare our results to theoretical prediction without losses; see S 8. A proper description of the long time dynamics will certainly require to account and understand such effects.

S5 Experimental uncertainties and inhomogeneities

Ideally, all atoms in the sample experience the same linear and quadratic Zeeman shift and the same quadratic light shift. However, in the experiment inhomogeneities from the magnetic field and light intensities lead to a spatial dependence of those shifts. An upper bound of the variation of Zeeman shifts can be deduced from RF-spectroscopy measurements done with bosonic erbium. From the width of the RF-resonance ($\approx 500\,$Hz) and the size of the cloud ($\approx 15\,\rm{\mu}$m) we estimate a maximum magnetic field gradient of $\leq 230\,$mG/cm, assuming the gradient as main broadening mechanism for the resonance width, neglecting magnetic field noise and Fourier broadening. This translates into a differential linear Zeeman shift of $\leq h\times 6\,$Hz between adjacent lattice sites in the horizontal x-y plane and $\leq h\times 12\,$Hz between adjacent planes along the z-direction. Together with the magnetic field values used in the spin dynamic experiments, the variation of the quadratic Zeeman shift is negligible compared to other inhomogeneities ($\leq h\times 0.1\,$Hz). The inhomogeneity of the quadratic light shifts can be estimated by considering the shape of the lattice light beams (Gaussian beams with waists of about $(w_{x},w_{y},w_{z})=(160,160,300)\,\rm{\mu}$m) and the resonance condition of the magnetic field, translated into a quadratic Zeeman shift of $h\times 71(1)\,\rm Hz$. This considerations can be used to obtain an estimation for a site dependent light shift compared to the center of the atomic sample. If we take a possible displacement of the atoms by $\leq 10\rm\mu m$ in all directions, from the center of the lattice to the center of the beams, into account, we can estimate an upper bound for the light shift of $\delta_{i}^{T}\leq h\times 2\,$Hz at $20$ lattice sites away from the center along the $y$ direction.

S6 Spin Hamiltonian

The experiment operates in a deep lattice regime, where tunneling is suppressed. At the achieved initial conditions, the 167Er atoms are restricted to occupy the lowest lattice band, and Fermi statistics prevents more than one atom per lattice site. In the presence of a magnetic field strong enough to generate Zeeman splittings larger than nearest-neighbor dipolar interactions, only those processes that conserve the total magnetization are energetically allowed de Paz et al. (2016). Under these considerations, the dynamics is described by the following secular Hamiltonian:

[TABLE]

Here the operators $F_{i}^{z,\pm}$ are spin $19/2$ angular momentum operators acting on lattice site $i$. The first two terms account for the site-dependent quadratic and linear shifts respectively, where $\delta_{i}$ includes both Zeeman terms and tensorial light shifts as discussed in the main text. $B_{i}=B+\Delta B_{i}$ denotes the linear Zeeman shift at site $i$. While the constant and uniform contribution, $B$, commutes with all other terms, thus can be rotated out, the spatially varying contribution, $\Delta B_{i}$, is relatively small in the experiment but still is accounted for in the theory calculations. The last term is the long-range dipolar interaction between atoms in different sites, with

[TABLE]

where $\theta_{ij}$ is the angle between the dipolar orientation set by an external magnetic field and the inter-particle spacing ${\bf r}_{i}-{\bf r}_{j}$. $V_{dd}\,\approx\,\frac{\mu_{0}g_{F}^{2}\mu_{B}^{2}}{4\pi d_{y}^{3}}$ corresponds to the interaction strength between two atoms, $i$ and $j$, separated by the smallest lattice constant $|{\bf r}_{i}-{\bf r}_{j}|=d_{y}=266$ nm and forming an angle $\theta_{{i},{j}}=\pi/2$ with the quantization axis. Here $g_{F}\approx 0.735$ is the Lande g-factor for Er atoms, $\mu_{0}$ is the magnetic permeability of vacuum and $\mu_{B}$ is the Bohr magneton. We compute $V_{dd}$ from

[TABLE]

where $\phi_{i}({r})$ denotes the lowest band Wannier function centered at lattice site ${i}$. For the experimental lattice depths $(s_{x},s_{y},s_{z})=(20,20,80)$ in units of the corresponding recoil energies, $V_{dd}$ is estimated to be $h\times 0.336\,$Hz.

S7 The GDTWA method

To account for quantum many-body effects during the dynamics generated by long-range dipolar interactions in these complex macroscopic spin $F=19/2$ 3D lattice array, we apply the so called Generalize Discrete Truncated Wigner Approximation (GDTWA) first introduced in Ref. Lepoutre et al. (2018). The underlying idea of the method is to supplement the mean field dynamics of a spin $F$ system with appropriate sampling over the initial conditions in order to quantitatively account for the build up of quantum correlations. For a spin-F atom $i$ with $\mathcal{N}=2F+1$ spin states, its density matrix $\hat{\rho}_{i}$ consists of $\mathcal{D}=\mathcal{N}\times\mathcal{N}$ elements. Correspondingly, we can define $\mathcal{D}$ Hermitian operators, $\Lambda^{i}_{\mu}$, with $\mu=1,...D$, using the generalized Gell-Mann matrices (GGM) and the identity matrix Bertlmann and Krammer (2008):

[TABLE]

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

[TABLE]

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

[TABLE]

for $1\leq\alpha<\mathcal{N}$

[TABLE]

With these operators, the local density matrix $\hat{\rho}_{i}$, as well as any operator $\hat{O}^{i}$ of local observables can be represented as

[TABLE]

and $\mu=1,2,...\mathcal{D}$. This allows expressing both one-body and two-body Hamiltonians in the form $\hat{H}_{i}=\sum_{\mu}c^{i}_{\mu}\Lambda^{i}_{\mu}$, and $\hat{H}_{ij}=\sum_{\mu,\nu}c^{ij}_{\mu\nu}\Lambda^{i}_{\mu}\Lambda^{j}_{\nu}$. The Heisenberg equations of motion for $\Lambda^{i}_{\mu}$ can be written as

[TABLE]

In the experiment, the initial state is a product state of single atom density matrices, $\hat{\rho}(t=0)=\prod\hat{\rho}^{i}(t=0)$. If we adopt a factorization $\langle\Lambda^{i}_{\mu}\Lambda^{j}_{\nu}...\Lambda^{k}_{\sigma}\rangle=\langle\Lambda^{i}_{\mu}\rangle\langle\Lambda^{j}_{\nu}\rangle...\langle\Lambda^{k}_{\sigma}\rangle$ for any non-equal $i,j,...k$ (i. e. each operator acts on a different atom) and arbitrary $\mu,\nu,\sigma$, Eq. S10 becomes a closed set of nonlinear equations for $\lambda^{i}_{\mu}=\langle\Lambda^{i}_{\mu}\rangle$. Within a mean-field treatment, the initial condition is fixed by $\lambda^{i}_{\mu}(t=0)={\rm Tr}[\Lambda^{i}_{\mu}\hat{\rho}(t=0)]$, which determines the ensuing dynamics from Eq. S10. This treatment neglects any correlations between atoms. In the GDTWA method, the initial value of $\lambda^{i}_{\mu}$ is instead sampled from a probability distribution in phase space, with statistical average $\overline{\lambda^{i}_{\mu}(0)}={\rm Tr}[\Lambda^{i}_{\mu}\hat{\rho}(t=0)]$. Specifically, each $\Lambda^{i}_{\mu}$ can be decomposed via its eigenvalues and eigenvectors as $\Lambda^{i}_{\mu}=\sum_{a^{i}_{\mu}}a^{i}_{\mu}\ket{a^{i}_{\mu}}\bra{a^{i}_{\mu}}$. We take $a^{i}_{\mu}$ as the allowed values of $\Lambda^{i}_{\mu}$ in phase space, then for an initial state $\hat{\rho}^{i}(t=0)$, the probability distribution is $p(a^{i}_{\mu})={\rm Tr}[\hat{\rho}^{i}(t=0)\ket{a^{i}_{\mu}}\bra{a^{i}_{\mu}}]$. From Eq. S10, each sampled initial configuration for the $N$ atom array, $\{a_{\mu}\}=\{a^{i1}_{\mu_{1}},a^{i2}_{\mu_{2}},...a^{iN}_{\mu_{N}}\}$ leads to a trajectory of $\Lambda^{i}_{\mu}$, which we denote as $\lambda^{i}_{\mu,\{a_{\mu}\}}(t)$. The quantum dynamics can be obtained by averaging over sufficient number of trajectories

[TABLE]

This approach has been shown capable of capturing the buildup of quantum correlations Schachenmayer et al. (2015b); Lepoutre et al. (2018).

S8 Incorporating experimental conditions in numerical simulation

In our experiment, the lattice filling fraction is not unity when the spin dynamics takes place. The reduced filling fraction is due to two effects: the finite temperature and atom loss during the initial state preparation. To account for the effect of a finite temperature, we first obtain the density distribution before ramping up the lattice from a Fermi-Dirac distribution $n^{0}({\bf r}_{i})=\frac{1}{1+{\rm exp}(\beta(\epsilon({\bf r}_{i})-\mu))}$, with parameters $\beta=1/k_{B}T$ and $\mu$ matching the inferred experiment temperature $T$, and the total atom number $N_{0}=2.4\times 10^{4}$. The function $\epsilon({\bf r}_{i})$ accounts for the weak external harmonic confinement. We compute the density distribution function after loading the atoms in the lattice, $n^{F}({\bf r}_{i})$ by simulating the lattice ramp which is possible since to an excellent approximation we can treat the system as non-interacting. Indeed, we neglect the dipolar interaction in the loading given that their magnitude is much lower than the Fermi Energy of the gas. In the numerical simulation, we then sample the position of atoms ${{\bf r}_{i}}$ in the lattice according to a distribution $p({\bf r}_{i})=n^{F}({\bf r}_{i})/N_{0}$. In practice, to reduce computation cost we need to reduce the total atom number in our calculations and use a smaller lattice with fewer populated lattice sites. In this case, we reduce the number of lattice sites by a factor $\xi=(N_{sim}/N_{exp})^{1/3}$, where $N_{sim(exp)}$ are the number of atoms in the simulation (experiment), while keeping the lattice spacings the same as in experiment, $(d_{x},d_{y},d_{z})=(272,266,544)\,$nm. That is, for an initial lattice with $L_{x}$ sites along $x$ direction, in our simulations there are $\xi L_{x}$ sites while the separation between two adjacent lattice sites is still $d_{x}$. We then sample the initial distribution of atoms in the lattice with $\tilde{p}(\tilde{\bf r}_{i})=\xi^{3}p(\xi\tilde{\bf r}_{i})$, which preserves the local density and is similar to sampling in a coarse-grained lattice. In our simulations, we chose $N_{sim}\gtrsim 350$ and checked that the convergence in $N_{sim}$ has been reached.

As discussed in Sec. S2, a fraction of atoms is lost during the ramp up and down of the magnetic field before initializing the spin dynamics over the sample. While a rigorous treatment on how these losses modify the distribution is not currently accessible with our current experimental setup, we try to account for it in the simulation by preferentially removing those atoms with a probablity $\propto p({\bf r}_{i})N_{\rm nn}$, where $N_{\rm nn}$ is the number of nearest neighbors (separation $\leq d_{y}$), until $N=\nu(0)N_{0}$ atoms are left. According to experiment estimates, the filling fractions before the initialization of the spin dynamics are $\nu(0)=0.6\sim 0.7$ (see Tab. S1 and S2). Figure S3 shows the histogram of neighbors in the resulting atom distribution. Such distribution effectively reduces the nearest-neighbor interactions and is found to give a better agreement with experiment.

Both the quadratic and linear shifts in the experiment are inhomogeneous across the lattice as discussed in Sec. S5, and we include them in our numerical simulation as site-dependent terms $\delta_{i}(\hat{F}_{i}^{z})^{2}$ and $B_{i}\hat{F}_{i}^{z}$, with $\delta_{i}=a|{\bf r}_{i}|^{2}$ and $B_{i}=b(x_{i}+y_{i}+z_{i})$. Based on experimental estimation, we have chosen the values of $a$ and $b$ such that $\delta_{i}=h\times 1.6\,$Hz ($h\times 0.7\,$Hz) at 20 sites along $y$ away from the lattice center, and $B_{i}$ differs by $h\times 6\,$Hz ($h\times 1.8\,$Hz) between adjacent sites, for Fig. 2  and 3 (Fig. 4) in the simulation.

S9 Short-time population dynamics

Considering a fixed initial atomic distribution over the lattice, the population dynamics at early times can be derived via a perturbative short time expansion

[TABLE]

Here the average $\langle\cdot\rangle$ is over the initial state, which is assumed to be a pure state, $\hat{n}_{m_{F}}=(\sum_{i}\mathcal{P}^{m_{F}}_{i})/N$, where $\mathcal{P}^{m_{F}}_{i}=|m_{F}\rangle_{i}{}_{i}\langle m_{F}|$ is the onsite projector for an atom at site $i$ in state $|m_{F}\rangle$ and $N$ denotes the total number of atoms. Note that here the sums are always carried out over the populated lattice sites in the initial lattice configuration. We obtain

[TABLE]

with

[TABLE]

[TABLE]

where $n_{m_{F}^{0}}$ denotes the population on the selected target state. To obtain Eq. S13, we have assumed that initially most of the population is in this target state, i.e. $n_{m_{F}^{0}}(0)\sim 1$. In the experiment, this assumption is always satisfied and therefore Eq. S15 is expected to reproduce well the short time dynamics.

The dependence of $\gamma(m_{F}^{0})$ on the initial state $m_{F}^{0}$ is a consequence of the dependence of dipolar exchange processes on the spin coherences, i. e.  $|\langle i:m_{F}^{0}+1,j:m_{F}^{0}-1|\hat{F}_{i}^{+}\hat{F}_{j}^{-}|i:m_{F}^{0},j:m_{F}^{0}\rangle|$. Therefore the smaller the value $|m_{F}^{0}|$ of the initial populated states, the faster the early time dynamics. Notably, up to order $t^{2}$ the initial dynamics is independent of quadratic shifts and external magnetic field gradients. This is because both of their corresponding Hamiltonians commute with the spin population operator $\hat{n}_{m_{F}}$. From this simple perturbative treatment one learns that by preparing different initial states with different $m_{F}^{0}$, the decay rates of the short time population dynamics provide information of $V_{\rm eff}$ and thus of the underlying dipolar couplings. As discussed in Sec. S2 and S8, the lattice filling fraction is not unity and the initial atomic density distribution in the lattice may vary from shot to shot. To account for this effect, we perform a statistical average of Eq. S14 calculated for each lattice configuration generated with the procedure in Sec. S8 to obtain the theoretical values in Fig. 3G and Fig. 4B.

It is important here to compare the predictions obtained from a simple mean-field analysis. In contrast to Eq. S15, neglecting quantum correlations yields

[TABLE]

At the mean field level therefore if initially the atoms are prepared such that $n_{m_{F}^{0}}(0)=1$, then there is no population dynamics. This is in stark contrast to the quantum systems where dynamics is enabled by quantum fluctuations. To extract $V_{\rm{eff}}$ from our experimental data and to compare it to the theoretical simulations we fit the initial dynamics with Eq. S15. We define the time scale for the fitting via $t_{\rm{fit}}<0.5\frac{\hbar}{V_{\text{eff}}}$, which corresponds to the timescale on which each atom did on average half a spin flip. We note that on this timescale the time evolution starts already to deviate from the short time expansion (Eq. 2), leading to a systematic downshift of the experimentally fitted $V_{\text{eff}}$; see Fig. 4B. However, a minimum timescale has to be chosen to ensure that the fit is performed using a large enough number of datapoints. Figure S4 shows exemplary the fit to the experimental data for $|m_{F}^{0}\rangle=|\text{--}9/2\rangle$.