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

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.

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 outof-equilibrium, with broad applications ranging from precision sensing and navigation, to quantum simulation and quantum information processing [1,2]. 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 [3,4]. 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 [5][6][7][8], trapped ions [9][10][11], polar molecules [12,13], and spin-3 bosonic chromium atoms [14][15][16], 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 nonequilibrium 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 longrange XXZ Heisenberg model. The simulator roots on the special atomic properties of 167 Er, 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 [17]. 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 [18], 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 [19]. 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 longrange 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 magnetizationconserving spin dynamics of our system can be conveniently written using spin-19/2 dimensionless angular momentum operatorsF i = {F x i ,F y i ,F z i }, acting on site i and satisfying the commutation relation [F x i ,F y i ] = iF z i . We use the eigenbasis ofF z denoted as |m F with 0 ≤ |m F | ≤ F [19][20][21]: is initialized by preparing all atoms in one starting state, here |-17/2 . We activate the spin dynamics by changing the magnetic field to set δ = 0. On early time scale dynamics are happening mainly among nearest neighbor atoms. Subsequently interactions between atoms at larger distances are involved in the dynamics.
The coupling constants V i,j = V dd d 3 y 1−3 cos 2 (θi,j ) r 3 ij , 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 r i and r j depends on their relative distance r ij = | r i − r j | and on their orientation, described by the angle θ i,j between the dipolar axis, set by the external magnetic field, and the interparticle axis; see Fig. 1B j terms describe dipolar exchange processes. The second sum denotes the single particle quadratic term δ i (F z i ) 2 with δ i = δ Z i + δ T i , accounting for the quadratic Zeeman effect ∝ δ Z i and tensorial light shifts ∝ δ T i . 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 [19]. 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 ≈ 10 4 Er atoms into a deep 3D optical lattice, following the scheme of Ref. [22]. 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 ∼ 10 Hz inside the planes and ∼ mHz between them [19]. The external magnetic field orientation, setting the quantization axis as well as the dipolar coupling strengths, is defined by the polar angles Θ and φ 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 [16].
Our experimental sequence to study the spin dynamics is illustrated in Fig. 1C. In particular, we prepare the system into the targeted m 0 F state by using the lattice-protection protocol demonstrated in Ref. [22]; see also Ref. [19]. At the end of the preparation, the majority of atoms are in the desired m 0 F (> 80%) at B ≈ 4 G. We note that atom losses during the spin preparation stage reduces the filling factor to about 60% of the initial one [19]. We then activate the spin dynamics by quenching the magnetic field to a value for whichδ = i δ 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 [19] 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 |-13/2 in Fig To benchmark our quantum simulator, we use a semiclassical phase-space sampling method, the so called generalized discrete truncated Wigner approximation (GDTWA) [16,19,[23][24][25]. 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 [19]. 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 [19]. To emphasize the beyond nearestneighbor 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 0 F 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 0 F , 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 ofF + iF − j , which gives a prefactor γ(m 0 . The initial dynamics can be well described by a perturbative expansion up to the second order [19], resulting in the analytic expression for the normalized population n m F (t) of the initial state: Here, ij is the overall effective interaction strength summed over N atoms and n m 0 F (0) denotes the purity of the initial state preparation. For a quantitative F from a fit to the experimental data (cyan triangles) and numerically computed from the initial spin distribution (black circles). Errorbars denote the 68% confidence interval of the fits. (H) All datasets with m 0 F < 0 used in (G) together with the corresponding full-GDTWA results (solid lines) plotted in units of the rescaled time τ = Veff/ · t. Note that all experiment and theory data are plotted for times, t ≤ 100 ms, of (A-F). To account for the different preparation fidelity, the populations of the initial states are shifted to 1 by adding a constant offset. For clarity error bars are omitted here.
analysis of the early-time spin evolution, we compare the theoretically calculated V 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 Veff estimated using the theoretically calculated V eff [26]. Figure 3G plots both, the theoretical and experimental V eff as a function of the initial m 0 F and highlights once more their quantitative agreement. The interaction parameter V eff , can also be used to rescale the time axis. As shown in Fig. 3H, all data sets now collapse onto each other for tVeff < 0.5, revealing the invariant character of the shorttime 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 Θ; see Fig. 1B. As exemplary shown in Fig. 4A for |-17/2 , the observed evolution speed of the spin populations strongly depends on Θ, changing by about a factor of 2 between Θ = 40 • and 80 • . The GDTWA results show a very good quantitative agreement with the experiment. We repeat the above measurements for different values of Θ and we extract V 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 (Θ c ≈ 35 • , φ c = 45 • ). We compare our measured V 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 Θ 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 [19].
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 [27] while adding Raman transitions to couple them, one could engineer nontrivial 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 [28][29][30].
We thank J. Schachenmayer for fruitful discussions and Arghavan Safavi-Naini for helping us understanding the loading into the lattice. We thank J. H. Becher and G. Natale for their help in the experimental measurements and for fruitful discussions. We also thank Rahul Nandkishore and Itamar Kimchi for reviewing the manuscript. The Innsbruck group is supported through an ERC Consolidator Grant (RARE, no. 681432) and a FET Proactive project (RySQ, no. 640378) of the EU H2020, and by a Forschergruppe  (1)) Hz, where ν ⊥ (ν ) are the trap frequencies perpendicular to (along) the horizontal ODT and ν 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 • 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 nm, d y = 266 nm, and d z = 544 nm. The lattice geometry is similar to the one used in our previous works [22,32]. 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 [22]. 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 ∈ x, y, z is given in the respective recoil energy E R,i with E R;x,y = h × 4.2 kHz and E R;z = h × 1.05 kHz, the atoms can be considered pinned on single lattice sites with low tunneling rates J x,y = h×10.5 Hz and J z = h × 1 mHz.

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 rapidadiabatic 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 [22]. 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 ms to 40.6 G to get a large enough differential quadratic Zeeman shift, which is on the order of ≈ h × 40 kHz. After the magnetic field ramp we wait for 80 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 |-19/2 into the |-7/2 state, we apply a radio-frequency (RF) pulse at 41.31 MHz and reduce the magnetic field with a linear ramp, e. g. by 500 mG in 42 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 ms to B = 3.99 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 × 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 |-19/2 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 |-17/2 , |-15/2 , |-9/2 , and |9/2 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 m F state to be 0.86(1) − 0.008(1) × 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.

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 [18] 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 θ p between the magnetic field axis and the axis of polarization of the light. Here, α t refers to the tensorial polarizability coefficient and ω to the light frequency. After the preparation of the respective spin state we start all our measurements at B = 3.99 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 Θ = 0 • for the different initial spin state and for the sets of measurements where |m F = |-17/2 and Θ ∈ (0 • , 80 • ). 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 nm and power of 7 W. This additional light is necessary because changing Θ reduces simultaneously θ p resulting in a smaller tensorial light shift and therefore in a shift of the resonance position to lower magnetic field values. For large Θ 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 Θ = 0 • but vary the initial m 0 F state we quench the magnetic field directly after the preparation, from 3.99 G to resonance. In contrast we use a different approach for the measurements where Θ is varied. After the preparation we ramp in 10 ms the additional laser beam to 7 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 x and B 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 ms and 1.2 ms, respectively. After holding on resonance for a certain time we quench the magnetic field back to 3.99 G and we rotate the latter back to Θ = 0 • . After a waiting time of 50 ms we perform a band-mapping measurement combined with a Stern-Gerlach technique, i.e. we ramp the lattice down in 1 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 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 , of its four neighbors, and of |-19/2 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 |-9/2 , |3/2 , and |5/2 . In the case of |3/2 residual atoms in |-19/2 , |-17/2 , and |5/2 are visible due to a non perfect preparation.

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 popula-tion 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 eff , t 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 0 F 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 ν = 0.58 for the |m 0 F = |1/2 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 0 F state, resulting larger for the central |m F 0 states. Similar numbers are obtained for the sets of measurement where we vary Θ (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.

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 (≈ 500 Hz) and the size of the cloud (≈ 15 µm) we estimate a maximum magnetic field gradient of ≤ 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 ≤ h × 6 Hz between adjacent lattice sites in the horizontal x-y plane and ≤ h × 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 (≤ h × 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) µm) and the resonance condition of the magnetic field, translated into a quadratic Zeeman shift of h × 71(1) 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 ≤ 10µm in all direc-tions, 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 δ T i ≤ h × 2 Hz at 20 lattice sites away from the center along the y direction.

SPIN HAMILTONIAN
The experiment operates in a deep lattice regime, where tunneling is suppressed. At the achieved initial conditions, the 167 Er 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 [15]. Under these considerations, the dynamics is described by the following secular Hamiltonian: Here the operators F z,± i 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 δ i includes both Zeeman terms and tensorial light shifts as discussed in the main text. B i = B + ∆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, ∆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 where θ ij is the angle between the dipolar orientation set by an external magnetic field and the inter-particle spacing r i − r j .
corresponds to the interaction strength between two atoms, i and j, separated by the smallest lattice constant |r i − r j | = d y = 266 nm and forming an angle θ i,j = π/2 with the quantization axis. Here g F ≈ 0.735 is the Lande g-factor for Er atoms, µ 0 is the magnetic permeability of vacuum and µ B is the Bohr magneton. We compute V dd from where φ 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 re-

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. [16]. 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 N = 2F + 1 spin states, its density matrixρ i consists of D = N × N elements. Correspondingly, we can define D Hermitian operators, Λ i µ , with µ = 1, ...D, using the generalized Gell-Mann matrices (GGM) and the identity matrix [33]: for α > β, 1 ≤ α,β ≤ N , With these operators, the local density matrixρ i , as well as any operatorÔ i of local observables can be represented aŝ and µ = 1, 2, ...D. This allows expressing both one-body and two-body Hamiltonians in the formĤ i = µ c i µ Λ i µ , and H ij = µ,ν c ij µν Λ i µ Λ j ν . The Heisenberg equations of motion for Λ i µ can be written as In the experiment, the initial state is a product state of single atom density matrices,ρ(t = 0) = ρ i (t = 0). If we adopt a factorization Λ i µ Λ j ν ...Λ k σ = Λ i µ Λ j ν ... Λ k σ for any nonequal i, j, ...k (i. e. each operator acts on a different atom) and arbitrary µ, ν, σ, Eq. S10 becomes a closed set of nonlinear equations for λ i µ = Λ i µ . Within a mean-field treatment, the initial condition is fixed by λ i µ (t = 0) = Tr[Λ i µρ (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 λ i µ is instead sampled from a probability distribution in phase space, with statistical average λ i µ (0) = Tr[Λ i µρ (t = 0)]. Specifically, each Λ i µ can be decomposed via its eigenvalues and eigenvectors as We take a i µ as the allowed values of Λ i µ in phase space, then for an initial stateρ i (t = 0), the probability distribution is p(a i µ ) = Tr[ρ i (t = 0) |a i µ a i µ |]. From Eq. S10, each sampled initial configuration for the N atom array, {a µ } = {a i1 µ1 , a i2 µ2 , ...a iN µ N } leads to a trajectory of Λ i µ , which we denote as λ i µ,{aµ} (t). The quantum dynamics can be obtained by averaging over sufficient number of trajectories This approach has been shown capable of capturing the buildup of quantum correlations [16,24].

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 (r i ) = 1 1+exp(β( (ri)−µ)) , with parameters β = 1/k B T and µ matching the inferred experiment temperature T , and the total atom number N 0 = 2.4 × 10 4 . The function (r i ) accounts for the weak external harmonic confinement. We compute the density distribution function after loading the atoms in the lattice, n F (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 r i in the lattice according to a distribution p(r i ) = n F (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 ξ = (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 exper- iment, (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 ξ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 withp(r i ) = ξ 3 p(ξr i ), which preserves the local density and is similar to sampling in a coarse-grained lattice. In our simulations, we chose N sim 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 ∝ p(r i )N nn , where N nn is the number of nearest neighbors (separation ≤ d y ), until N = ν(0)N 0 atoms are left. According to experiment estimates, the filling fractions before the initialization of the spin dynamics are ν(0) = 0.6 ∼ 0.7 (see Tab. S1 and S2). Figure S3 shows the histogram of neighbors in the resulting atom distribution. Such distribution effectively reduces the nearestneighbor 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 sitedependent terms δ i (F z i ) 2 and B iF z i , with δ i = a|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 δ i = h × 1.6 Hz (h × 0.7 Hz) at 20 sites along y away from the lattice center, and B i differs by h×6 Hz (h×1.8 Hz) between adjacent sites, for Fig. 2 and 3 (Fig. 4) in the simulation.

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 (S12) Here the average · is over the initial state, which is assumed to be a pure state,n m F = ( i P m F i )/N , where P m F i = |m F ii m F | is the onsite projector for an atom at site i in state |m F 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 with where n m 0 F 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 0 F (0) ∼ 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 γ(m 0 F ) on the initial state m 0 F is a consequence of the dependence of dipolar exchange processes on the spin coherences, i. e. | i : m 0 F + 1, j : m 0 F − 1|F + iF − j |i : m 0 F , j : m 0 F |. Therefore the smaller the value |m 0 F | 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 operatorn m F . From this simple perturbative treatment one learns that by preparing different initial states with different m 0 F , the decay rates of the short time population dynamics provide information of V 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, ne- At the mean field level therefore if initially the atoms are prepared such that n m 0 F (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 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 fit < 0.5 Veff , 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 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 0 F = |-9/2 .